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CHAPTER  I 


INTRODUCTION 

The  object  of  this  work  is  to  use  numerical  methods  to 
solve  electromagnetic  field  problems  involving  scattering 
from  rather  arbitrarily  shaped  inhomogeneous  penetrable 
bodies.  While  our  aim  is  to  analyze  a missile  in  the  pre- 
sence of  an  electrically  inhomogeneous  exhaust  plume,  the 
techniques  discussed  here  are  useful  in  other  areas  of  elec- 
tromagnetics, such  as  scattering  by  rain  drops,  power 
absorption  in  biological  tissues,  dielectric  lenses,  etc.  A 
primary  requirement  of  any  numerical  method  is  that  the 
technique  should  be  capable  of  simulating  the  actual  physical 
situation  as  closely  as  possible,  while  simultaneously  pro- 
viding an  efficient  method  of  solution.  For  a numerical 
study  of  scattering  by  inhomogeneous  dielectric  bodies,  one 
must  choose  among  a variety  of  techniques,  all  of  which  can 
be  said  to  fall  into  one  of  the  following  two  categories:  (a) 
Integral  equation  formulations  and  (b)  Differential  equation 
methods.  The  usual  surface  and  volume  integral  equation 
formulations  with  numerical  solution  by  the  method  of  moments 
[1]  and  the  extended  boundary  condition  approach  [2]  fall 


into  the  category  of  (a)  while  finite  difference  methods  and 
the  unimoment  method  [3]  fall  into  (b).  Thus  a wide  variety 
of  approaches  are  available.  Some  of  the  main  features  of 
each  of  these  methods  are  given  in  the  following. 

The  volume  integral  equation  is  based  on  relating  the 
polarization  current  in  terms  of  the  total  field,  comprising 
the  incident  field  and  the  scattered  field.  By  associating 
an  unknown  current  coefficient  with  each  point  inside  the 
region,  the  integral  equation  is  converted  into  a matrix 
equation  which  can  be  easily  solved  for  the  unknown  coef- 
ficients. Since  the  region  of  the  scatterer  is  represented 
point  by  point,  an  arbitrary  inhomogeneity  and  shape  is 
easily  handled  in  this  approach.  The  approach,  however, 
leads  to  very  le-Ae  matrices  which  makes  the  method  unat- 
tractive due  to  the  limited  core  storage  on  the  computer. 

The  integral  equation  approach  is  well  suited  either  for 
homogeneous  penetrable  bodies  or  for  a body  either  modeled 
by  or  made  up  of  layers  of  homogeneous  regions.  The  usual 
procedure  in  this  case  is  to  set  up  the  coupled  integral 
equations  in  ter-s  of  equivalent  electric  and  magnetic  cur- 
rents on  the  surfaces  of  the  homogeneous  region.  By  expand- 
ing the  unknown  currents  in  terms  of  suitable  basis  functions 
and  adopting  suitable  testing  functions,  the  coupled  integral 
equations  are  reduced  to  a matrix  equation  for  the  unknown 
coefficients  of  the  basis  functions. 
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For  a body  made  up  of  a large  number  of  layers  or  for 
an  inhomogeneous  body  modeled  as  being  made  up  of  layers  of 
homogeneous  regions,  the  above  approach  can  lead  to  very 
large-sized  matrices  due  to  a simultaneous  solution  of  the 
equivalent  currents  on  all  the  layers.  Since  the  fields 
induced  in  any  region  between  two  layers  are  due  to  the 
equivalent  currents  on  adjacent  layers,  the  resulting  matrix 
is  block  tridiagonal.  This  property,  as  recognized  by 
Pogorzelski  (4],  yields  an  iterative  procedure  for  solving  for 
the  currents  on  the  outermost  layer  in  terms  of  the  currents 
on  the  inner  layers.  Such  an  iterative  procedure  has  the 
advantage  that  the  sizes  of  the  matrices  involved  in  the 
iteration  are  much  smaller  than  the  overall  matrix  size 
that  would  be  required  if  the  currents  on  all  layers  were  to 
be  solved  simultaneously.  The  surface  equivalence  approach 
is,  in  principle,  applicable  to  all  (layered)  inhomogeneous 
scatterers,  regardless  of  shape. 

The  extended  boundary  condition  approach  proposed  by 
Waterman  (2]  expresses  the  fields  in  terms  of  integrals  over 
surfaces  separating  the  homogeneous  regions  around  the  scat- 
terer.  However,  one  uses  here  the  fact  that  in  all  regions 
complementary  to  those  in  which  the  equivalence  is  valid,  the 
fields  must  vanith.  Within  these  null  field  regions,  the 
integral  expressions  for  the  fields  are  expanded  in  spherical 
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fey lindrical)  harmonics  of  the  standing  wave  type  for  in- 
terior three-dimensional  (two-dimensional)  regions  and  of 
the  outgoing  type  for  exterior  regions.  Since  the  fields 
vanish,  the  coefficients  of  the  harmonics  must  also  vanish 
and  one  obtains  a set  of  equations,  each  of  which  involves 
integrals  over  the  equivalent  sources  on  the  surface.  The 
surface  sources,  both  electric  and  magnetic,  are  usually  ex- 
panded in  spherical  (cylindrical)  harmonics  as  well  and  this 
eliminates  one  of  the  surface  sources,  thus  reducing  the 
number  of  unknowns.  In  the  literature,  the  extended  boundary 
condition  approach  is  also  known  as  the  T-matrix  approach. 
Peterson  and  Strom  [5]  have  extended  this  approach  to  multi- 
layered dielectric  scatterers  and  Strom  [6]  has  further 
extended  the  approach  to  multiple  inhomogeneous  scatterers. 
Since  the  method  depends  on  the  object  having  an  interior 
region  in  which  a circumr cr ibed  sphere  (circle)  can  be 
placed,  the  method  is  better  suited  for  nearly  spherical 
(cylindrical)  bodies  chan  for  thin  scatterers. 

The  unimoment  technique  developed  by  Mei  [3]  and  ex- 
tended by  Chang  [7]  and  Morgan  [8]  essentially  studies  the 
scattering  problem  through  a differential  equation  formula- 
tion. According  to  this  approach,  a spherical  (cylindrical) 
region  surrounds  the  three-dimensional  (two-dimensional) 
scatterer.  The  minimum  radius  of  this  region  should  be  so 
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as  to  totally  enclose  the  scatterer  so  that  the  scattered 
fields  in  the  exterior  region  can  be  expanded  in  terns  of 
outgoing  spherical  (cylindrical)  harmonics.  The  incident 
field,  of  course,  may  be  expanded  in  terms  of  incoming  har- 
monics. A wave  equation  for  an  appropriate  field  quantity 
is  next  solved  for  in  the  interior  region.  The  boundary 
conditions  for  the  tangential  fields  are  then  enforced  across 
the  spherical  (cylindrical)  boundary.  This  results  in  a set 
of  equations  which  determine  the  coefficients  of  the  unknown 
scattered  fields.  As  one  notes,  the  major  effort  Involved 
in  this  approach  is  in  solving  the  differential  equation  in 
the  interior  region.  Either  a finite  difference  approach  or 
the  finite  element  method  [9],  a numerical  approach  for  solv- 
ing differential  equations  that  has  been  highly  developed  by 
structural  engineers,  can  then  be  utilized  to  solve  this 
differential  equation.  In  the  latter  approach,  the  interior 
region  is  typically  divided  into  a number  of  triangular 
sections  called  elements.  Over  each  element  the  field  is 
represented  by  suitable  expansion  functions  that  express  the 
field  within  an  element  as  a function  which  interpolates  the 
value  of  the  field  at  the  nodes  of  the  element.  By  minimiz- 
ing a stationary  formula  associated  with  the  differential 
operator  with  respect  to  the  nodal  coefficients,  one  obtains 
a matrix  equation  for  the  nodal  field  values  on  the  interior 
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in  terms  of  the  field  values  at  the  nodes  lying  on  the 
spherical  (cylindrical)  boundary.  Then  using  the  spherical 
(cylindrical)  modes  as  basis  functions  to  excite  the  electric 
fields  at  the  boundary,  one  obtains  interior  fields  for  each 
distinct  mode  of  the  harmonics.  Orthogonality  of  the  basis 
functions  enables  one  to  determine  the  interior  fields  in 
terms  of  the  scattered  fields.  Equating  the  exterior  and 
interior  tangential  magnetic  fields  at  the  boundary  yields  a 
matrix  equation  which  determines  the  scattered  field  coef- 
ficients. An  inherent  feature  of  the  matrix  so  obtained  is 
that  it  is  banded  and  an  efficient  utilization  of  this  prop- 
erty allows  one  to  soive  problems  involving  a rather  large 
number  of  unknowns.  Furthermore,  just  as  one  can  express 
the  interior  fields  in  terms  of  interpolatory  functions  over 
each  element,  one  can  similarly  approximate  the  spatial 
variation  of  the  physical  parameters  viz . , the  permittivity 
and  permeability,  by  means  of  the  same  interpolatory  func- 
tions. This  latter  feature  of  the  unimoment  method  enables 
one  to  solve  for  fields  from  arbitrary  inhomogeneous  (i.e. 
not  necessarily  layered)  scatterers.  As  with  the  extended 
boundary  condition  approach,  the  unimoment  method  is  more 
suitable  for  scatterers  which  are  almost  circular  or  spher- 
ical in  shape. 

From  the  above  discussion  it  is  apparent  that  the  choice 
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of  a method  should  be  made  from  the  point  of  view  of  sim- 
plicity, accuracy,  and  efficiency  of  the  method  as  it  applies 
to  the  geometry  and  scatterer  in  question.  In  order  to  more 
clearly  define  some  of  the  considerations  involved  In  the 
choice  of  the  method,  we  consider  the  application  of  the  unimo- 
ment method  to  scattering  from  two-dimensional  layered  di- 
electric cylinders  in  Chapter  11.  We  also  note  some  of  the 
advantages  of  a local  as  compared  to  a global  coordinate 
formulation.  In  Chapter  II  the  unimoment  method  approach  is 
compared  with  the  iterative  solution  procedure  for  the  sur- 
face integral  formulation  for  scattering  by  layered  dielec- 
tric cylinders.  Chapter  III  deals  with  the  application  of 
the  iterative  solution  procedure  for  the  surface  integral 
equation  to  layered  bodies  of  revolution.  In  Chapter  IV,  the 
approach  is  extended  to  treat  missile-plume  problems,  which 
are  also  reduced  to  a block-tridiagona 1 form. 

In  the  course  of  this  work,  the  equivalence  principle 
is  used  extensively.  One  normally  uses  the  equivalence 
principle  to  set  up  coupled  integral  equations  for  unknown 
electric  and/or  magnetic  currents.  In  Appendix  A,  different 
types  of  integral  equation  formulations  are  considered.  The 
discussions  there  parallel  and  extend  slightly  the  work  in 
this  area  by  Harrington  and  Mautz  (lOj. 
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CHAPTER  II 

SCATTERING  BY  INHOMOGENEOUS  DIELECTRIC  CYLINDERS 

In  this  chapter,  we  present  a comparative  study  of  the 
unimoment  method  and  the  iterative  solution  procedure. 

Since  the  principal  features  of  the  methods  when  applied  to 
a general  problem  can  be  illustrated  through  specific  exam* 
pies  for  which  alternative  solutions  are  available,  we 
consider  herein  scattering  by  layered  dielectric  cylinders. 

In  Section  2.1  a brief  description  of  the  unimoment  method 
is  given.  While  we  follow  essentially  the  approach  of  [7], 
we  do,  however,  present  a simpler  representation  of  the 
unimoment  matrix  than  is  found  there.  Sec.  2.2  gives  some 
insight  into  the  iterative  solution  of  the  surface  integral 
method  which  leads  to  a b lock-tr idi agonal  moment  matrix.  A 
comparison  with  the  unimoment  method  is  then  made  to  point 
out  the  applicability  and  limitations  of  the  .wo  techniques 
for  general  problems  involving  inhomogeneous  dielectric 
bodies. 

2 . 1 Unimoment  Method  for  Scattering  from  Dielectric  Cylinders 
Fig.  2.1  shows  an  arbitrary  cylindrical  scatterer  upon 
which  a plane  wave  is  incident.  We  shall  restrict  our 


9 


Figure  2.1.  Geometry  of  the  scattering  problem. 
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discussions  herein  to  a TM  polarized  incident  wave.  The 
reader  is  referred  to  [7]  for  the  TE  polarization  case. 

The  scatterer  is  enclosed  by  a circle  of  radius  a.  As- 
suming that  the  axis  of  the  cylinder  is  in  the  z direction, 

let  E**(r,<}>)  be  the  total  field  inside  the  circle,  Einc(r,<j>) 
z z 

s c 

the  incident  field»and  E (r,<ji)  the  scattered  field  outside 
the  circle.  The  tangential  components  of  the  fields  must  be 
continuous  across  the  circular  boundary.  Thus,  we  have 


E^1  (a,  <j>)  - E*nC(a,<j>)  + E®C(a,<|>),  (2.1) 


»inc 

Kc 

3r  1 

3r 

r ■ a 

+ 3r 

r ■ a 

Since  Ez  (r,t}>)  is  a scattered  field,  we  express  it  in 
terms  of  outgoing  cylindrical  harmonics  as 

00 

E^C(r,4>)  * y*  r)i.\  cos  mji  + B sinn<j>[  , 

z n o n n J 

n ■ o (2.3) 

(2) 

where  H (x)  is  the  Hankel  function  of  the  second  kind, 
n 

k ■=  to  , is  the  free  space  propagation  constant, 

and  are  arbitrary  constants  to  be  determined.  Since  the 
incident  field  can  also  be  expanded  in  a series  of  cylin- 
drical harmonics,  let  us  therefore  express  the  field  of 


' »-=r*vr-i 
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region  II  evaluated  on  the  circular  boundary  in  terms  of 
cylindrical  harmonics.  Thus,  we  have 


.11/  .. 
c (a , <p, 

z 


(C  cos  n<{>  + D sin  n ) . (2.4) 

n n 


n * o 


Each  Fourier  harmonic  in  (2.4)  can  be  thought  of  as  repre- 


senting the  evaluation  on  the  boundary  of  one  term  of  a 

(C) 

complete  set  of  linear  independent  partial  fields,  (r,<{>) 

which  satisfy 


V2^g)( r,$)  + k2er(r,$)i|;^s)  (r,4>)  » 0, 


and  the  boundary  conditions 


(a , <{>) 


cos  n$  , 


(2.5) 


(2 . 6a) 


i ' 
i ! 


sin  n4>  . 


Thus  the  interior  field  is  given  by 


n n 


n n 


n = o 


.<§) 


(2.6b) 


E^1  (r, $)  - ^ + D_<J»*(r,<|>)  - (2.7) 


We  may  determine  (r,<p)  with  the  above  boundary  conditions 

by  solving  (2.5)  through  any  of  the  standard  techniques  of 


i ' 


r • 'VS”^^ — 


12 


solving  the  second-order  differential  equation.  We  shall  be 
adopting  the  finite  element  method,  a discussion  of  which  is 
reserved  until  the  next  section.  Using  (2.3)  and  (2.4)  in 
(2.1),  we  get 


2 


n - o 


(C  cos  n<j>  + D sin  n 4> ) - E*nc(a,<J>) 
n n z 


\*H(2)(k  a)  [A  cos  n$  + B sin  n<M  (2.8) 

y n o 1 n r n 


n ■ o 


I [‘ 


n - o 


the  use 

of  (2.3)  and  (2. 

7)  in  (2.2),  we 

obtain 

r 

. Einc 

n 

+ n 2. 

1 - 1 

C 

3 

OJ 

*-» 

“ u . 

n 9r 

9r 

I 

*■  l 

r - a 

r - a J 

r * a 

+ k H^2^  (k  a)  [A  cos  n<J>  + B sin  n^]  . (2.9) 

^ on  o 1 n n 


n ■ 1 


Invoking  standard  orthogonality  relationships  we  obtain  from 

(2.8) 


C 

n 


A H(2) (k  a)  + fc 
n n o n 


(2.10a) 


D 

n 


B H(2)(k  a)  + fs 
n n o n 


(2.10b) 


where 
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f 


c 

n 


f 


8 

n 


(a, <p)  cos  n<j>  d<J>, 


(a,if>)  sin  n<J>  d<}>  . 


(2.11) 


(2.12) 


Substituting  (2.10)  - (2.12)  in  (2.9),  we  obtain  finally 


I [A»  {<2><k-*>  "" 


o'  3r 


n * o 


r - a 


( 2)  1 

- k H'  (k  a)  cos  n<J) 
o n o 


+ »„  K2)<V>  IT  t.a  - ko“"),(V>  *1”  "♦}] 


3t’ 


,inc 


r « a 


+ f 


n M o 


r » a 


3ij>s 
s r n 

n 3r 


r * a ^ • 


(2.13) 


In  computations,  the  summation  in  (2.13)  must  be  truncated 
to  N terms.  The  choice  of  N is  generally  slightly  greater 
than  kQa.  To  obtain  a matrix  equation  for  the  unknown 
coefficients  An,  Bn>  (2.13)  may  be  multiplied  by  cos  n 4> 

(n  = 0,  1,  2...N)  and  sin  n<t»(n  = 1,  2...N),  integrated  over 
the  interval  0 to  27r.  (Note  that  this  is  equivalent  to 


<5' 
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expanding  the  continuity  equation  (2.2)  into  Fourier  modes). 
Once  A^,  have  been  computed  the  scattered  field  is  easily 
computed  from  (2.3).. 

2 . 2 Interior  Problem 

The  interior  field  satisfies  the  differential  equation 
(2.5).  Since  the  exact  solution  to  the  inhomogeneous  equa- 
tion is  possible  only  for  limited  types  of  inhomogeneity, 
we  resort  to  the  so-called  finite  element  method  [9]  for 
solving  for  In  this  approach,  the  solution  for  is 

obtained  by  minimizing  a variational  functional  associated 
with  the  differential  equation.  The  interior  of  the  circle 
of  radius  a is  usually  divided  into  a number  of  triangular 
subregions  which  are  known  as  1 elementc".  The  function 
is  expanded  over  rach  element  in  terms  of  suitable  functions 
called  "trial  functions".  The  values  of  the  trial  functions 
are  specified  at  certain  points  (nodes)  on  the  triangles. 
Typically  these  nodes  are  at  the  vertices  of  the  triangular 
elements,  but  in  higher  order  schemes  [9],  may  also  be  at, 
say,  the  mid-points  of  the  sides  of  the  triangles.  By  using 
the  trial  functions  over  each  element  and  minimizing  the 
functional  with  respect  to  the  nodal  values  within  an  ele- 
ment, one  obtains  a matrix  equation.  This  matrix  equation 
may  be  then  solved  to  determine  at  the  nodal  points. 
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For  TM  polarized  waves,  1 1>  satisfies 

n 

V2«Pn(r,<}»)  + k^er(r,4»)tjjn(r,<))  = 0, 

where  ifi  (r,<p)  » E (r,$).  The  solution  to  this  equation  is 
the  same  as  one  would  obtain  by  minimizing  the  functional 


1 


Ki2  ~ ^oer(r»^>^(r»<^> 


(2-14) 


where  S corresponds  to  the  region  over  which  (2.5)  is  valid. 

We  divide  the  region  S into  a number  of  triangluar  regions 

(Fig.  2.2).  Let  Y , i - 1,  2,  ...  K + K',  be  the  nodal 

i 

values,  where  K is  the  total  number  of  nodes  inside  the 
region  S and  K’  is  the  total  number  of  nodes  on  the  artificial 
circle.  Approximate  ^ by  a linear  function  over  each  tri- 
angle. Minimizing  the  functional  3.  in  (2.14)  with  respect 

to  Y , i - 1,  2...K,  one  obtains 
ni 

QU  - TU*  , (2.15) 


where 
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K+l 


K+2 


°K+K ' 


and  Q and  T are  the  coefficient  matrices  obtained  from  (2.14) 
over  each  element.  Since  ¥ , i - K+l,  K+2,  . ...K+K',  are 


known  for  each  mode  n,  one  can  solve  for  U from  (2.15).  Once 

3<f>, 

is  evaluated  along  the  artificial 

r - a 


to  is  known,  -r— 
n dr 


circle.  Having  obtained  the  numerical  derivative  of  to  on 

n 

the  circle  (r«a),  we  may  then  evaluate  the  scattered  field 

coefficients  A , B from  (2.13). 
n n 

The  remaining  problem  is  thus  the  evaluation  of  the 
elements  of  the  matrices  Q and  T.  It  is  in  this  step  that 
we  follow  a slightly  different  derivation  than  that  of  Chang 
and  Mei  [7).  If  <j>  - toR  is  the  field  in  any  element  cor- 
responding to  the  nC^  mode,  Chang  and  Mei  express  the 
assumed  linear  variation  of  this  field  in  terms  of  a fixed 
global  (x,y,z)  coordinate  system  as 


- ax  + by  + c , 


(2.16) 


where  a,  b,  and  c are  the  expansion  coefficients  to  be 
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determined.  If  (1,2,3)  are  the  nodes  of  the  triangular 
element  and  <j>^,  <f> ^ » 4>.j  are  the  respective  nodal,  values , then 
the  expansion  coefficients  a,  b,  c for  this  element  can  be 
easily  found  in  terms  of  ^ as 


-a 


a 

xi  *i  f 

V 

b 

m 

*2  >-2  1 

*2 

c 

.*3  >'3  X. 

♦3. 

(2.17) 


An  alternative  and  simpler  representation  for  <{>  is  in  terms 
of  a local  coordinate  system,  wherein  <J>  is  expressed  in  terms 
of  area  coordinates  {9].  Referring  to  Fig.  2.3,  the  field  in 
the  element  is  written  in  terms  of  area  coordinates  as 

<t>  « + <*>2A2  + ^3A3)’  (2.18) 

where  A is  the  total  area  of  the  triangle  and  A , m * 1 , 2 , 3 

m 

are  the  sub-areas  shown  in  Fig.  2.3.  It  should  be  noted 
herein  that  the  global  coordinate  and  local  coordinate  rep- 
resentations are  two  different  ways  of  expressing  the 
assumed  linear  variation.  The  two  representations  themselves 
can  be  related  to  one  another.  The  relationship  between  the 


two  can  be  obtained  as 


20 


Y 

‘l 

1 

1 " 

rv 

A 

x 

m 

X1 

X2 

x3 

a2 

T 

y 

» ■ 

yl 

* 

y 2 

y3 

a3 

A 

» 

HJ* 

1 

P 2 3 

Y23 

■ 

“X23 

1 

A, 

4 

A 

m 

P31 

Y31 

"X31 

X 

A3 

J 

A 

P12 

Y12 

“X12 

y 

m m 

wher'  p»»  ■ *i  - 2k  'v,  • v.>- 

T»n  * Yi  ■ H l!,»  - J’n1 


(2.19) 


(2.20) 


X 


an 


x^]  , £.,m,n  - 1,2,3 


and  (x^,y^)  , i * 1,  2,  3,  are  the  coordinates  of  the  nodes 
of  the  element.  Since  the  integrand  of  the  functional  I in 
(2.14)  involves  derivatives,  we  note 


3x 


[ 


f 


3 


i * 1 


1 

A 


- UL 

i 3A , 


, (2.21) 
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_9_ 

9y 


[^.^^)]  - i -ki^ 


1 = 1 

The  functional  I in  (2.14)  over  the  element  is 


123 


ffh2-  ko£r^2J  ds  • 


123 


We  note  that 


Vd>  « x +•  y 

vv  9x  x 9y  y 


Using  (2.18),  (2.21)  and  (2.22),  we  get 


A 

fi  ♦ 

V 1 Y 

L A Yi 

X T 

2,  A Xi 

U - 1 J 

i = 1 

Hence 


|v<M 


3 A 

CM 

('  * ) 

V iv 

X 

V J, 

2L  A Yi 

T 

Z * i 

i « 1 

[i  - i 

Minimizing  I.  with  respect  to  one  has 

123 


91 


123  9 


9 3 <f> 


II 


[ | V d>  | 2 - k^e^2]  ds  . 


‘123 


(2.22) 


(2.23) 


(2.24a) 


(2.24b) 


(2.23) 
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Substituting  (2.24)  and  (2.18)  into  the  above  and  noting 
the  following  integrals. 


■§  , i - 1,2,3 


i * j 


12  3 


one  obtains 


31, 


123 


3$! 


V *m 
y y 

Y,  + 2 

y •—  x 

^ A m 

1 

4m  A m 

■i 

a «*  1 

k^AC 

- <2*i + *2 + *3}  * 


(2.26) 


Herein  we  have  assumed  that  er  is  constant  over  the  element. 
If  varies  over  the  element,  then  one  can  expand  in 
terms  of  a suitable  polynominal  (in  the  area  coordinates) 
and  proceed  as  above  to  evaluate  the  second  term  in  (2.25). 
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If  there  are  p triangles  around  node  1,  then  one  has  to 
evaluate  the  integral  in  (2.25)  over  each  of  these  p tri- 
angles to'obtain'a'  single  element'of  the  matrix  Q’  or"  Tl 
From  the  above  it  is  clear  that  a proper  bookkeeping  of  the 
nodes  that  make  up  an  element  of  the  finite  element  mesh  is 
mandatory.  However,  it  is  the  nodes  which  correspond  to  the 
elements  of  the  matrices  Q and  T.  Accordingly,  it  is  the 
interconnecting  nodal  arrangement  which  plays  a key  role  as 
compared  to  the  nodes  that  make  up  each  element.  Computa- 
tionally, one  could  search  through  a matrix  which  lists  the 
nodes  making  up  an  element  and  find  which  elements  of  the 
mesh  contribute  to  which  elements  of  the  matrix.  This 
appears  to  be  the  procedure  adopted  in  [7],  An  alternative 
approach  that  is  computationally  more  efficient,  however,  is 
to  define  the  elements  through  a connection  matrix  N * [ n ^ ^ , 
whose  elements  are  the  nodal  numbers  of  the  j ^ node  (num- 
bered counter-clockwise) ■ connected  to  node  i (see  Fig.  2.4). 
Such  ar  approach  has  the  advantage  of  immediately  identifying 
the  interconnecting  nodes.  It  also  makes  the  numbering  of 
the  elements  superfluous,  as  we  shall  see.  Expanding  out 
the  terms  in  (2.26),  we  obtain 
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A. 


-HI  = ^ (x2  2 

3$.  A 'A1  + VV1 


+ | [ CXXX2  + Y1Y2)(j)2  + | (X1X3  + Y1Y3)<*>3] 


k2e 
o r 


3 [A*1  + f *2  + | *3]  . 


(2.27) 


Define  the  vector  location  of  node  i,  with  coordinates  (x^, 

yt)  as 


ri  • xi  + Jyt 


and  the  vector  from  node  i to  node  j as 


ij 


- r , 


Note  that 


r 


12 


x 


r23  * 


31 


r , x 

31 


12 


With  the  above  convention,  (2.27)  reduces  to 
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31 
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^2  + 
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12 


2A 


k2e 
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I A<J> ! 


(2.28) 


wherein  the  various  vectors  are  shown  in  Fig.  2.5  . We  note 
here  that  the  self-node  contribution  of  the  Laplacian  term 
is  proportional  to  the  ratio  of  the  square  of  the  distance 
of  the  side  opposite  to  a node  to  the  area  of  the  element, 
while  the  coefficients  of  the  mutual  terms  in  the  Laplacian 
are  proportional  to  the  cotangent  of  the  angle  formed  by 
the  sides  at  the  opposite  non-self  node 

One  can  obtain  the  Laplacian  terra  above  in  an  alterna- 
tive way.  Expanding  $ In  area  coordinates  as  in  (2.18)  and 

1 

noting  that  7*  is  given  by  (2.24a),  one  can  use  7 = V • V<)> 

avid  approximate  the  divergence  as 


V 


(2.29) 


where  n is  the  normal  to  the  contour  C (in  the  plane  of  C) 
enclosing  the  triangles.  Using  (2.24a)  in  (2.29),  one 
obtains  exactly  the  Laplacian  term  given  in  (2.27).  Such 
an  observation  is  very  useful  when  dealing  with  surfaces 
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that  may  be  irregularly  shaped  or  for  problems  wherein  the 
functional  I cannot  be  easily  formulated. 

In  general,  if  there  are nodes  around  node  i,  as  shown 
in  Fig.  2. A,  differentiating  the  ve-iational  form  with  re- 
spect to  each  of  the  nodal  values  yields 


i - 1,2. . .K,  (2.30) 


wher 


-i 

r.  » r - r , 
j n.  i 

J 


r.,,  , * r - r 

j+1.3  ni , j+i  "i 


J 


. 'SgjHpCBH 


A1  * — i • (r*  x r1-) 

J 2 1 J+l  V 


j = 1 + (j  - 1) 


mod  IT 


Equating  (2.30)  to  zero  so  as  to  satisfy  the  stationarity 

property,  we  obtain  the  elements  q and  t of  the  matrices 

nn>  nm 

Q and  T of  (2.15)  as 


j - 1 


l l%t±il ’ - l *;  • 

j - 1 


(2.31) 


-n  . rn  , ?n  . ?n 

a.m-1 a-1  + m+l.m  a 

2Ci  ^ 


- r~(An  + An  )i  * “ ' 

l2U»  + VlJ  ’ ’ m > K • ' 


(2.32) 


Using  (2.31)  and  (2.32),  we  may  calculate  the  interior  fields 
from  (2.15)  for  each  mode  specified  on  the  artificial  circle. 
Once  the  interior  fields  are  known,  a finite  difference 

a*n  dK 

scheme  yields  -r — . With  t|>  and  -r — known,  the 

dr  r * a n Dr  r =*  a 

scattered  field  coefficients  A , B are  determined  from 

n n 

(2.13) . 
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2 . 3 Numerical  Results  of  Unimoment  Method 

A computer  program  was  first  written  using  the  procedure 
In  [7]  which  uses  the  global  coordinate  representation.  Lat- 
er, a separate  program  using  the  local  coordinate  representa- 
tion discussed  here  was  written.  Accuracy  of  the  programs 
was  checked  by  comparing  with  the  exact  solutions  for  a di- 
electric circular  cylinder.  Figs.  2.6  and  2.7  show  the 
scattered  field  patterns  for  a two-layered  circular  cylinder 
and  a two-layered  elliptical  cylinder,  respectively.  The 
agreement  between  the  exact  solution  in  Fig.  2.6  and  the 
moment  method  solution  in  Fig.  2.7  is  quite  good.  The 
minor  difference  one  notices  in  the  backscatter  direction 
Is  attributed  primarily  to  inaccuracies  introduced  by  equa- 
ting the  analytically  exact  normal  derivative  of  the 
exterior  field  with  the  numerically  derived  value  for  the 
interior  field  which  was  computed  by  a backward  difference 
at  the  boundary. 

A few  additional  comments  concerning  the  use  of  the 
local  coordinate  representation  are  in  order.  The  use  of 
the  connection  matrix  simplifies  inputting  the  data  to  a 
computer  code  as  well  as  saves  significantly  the  time  spent 
in  searching  through  a list  of  element  numbers  to  determine 
the  connecting  nodes.  This  fact  becomes  more  apparent  if 
one  realizes  that  the  element  numbers  are  superfluous  as 


er 
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compared  to  the  nodes  that  make  up  the  elements.  Further, 
the  use  of  the  area  coordinate  representation  allows  one  to 
explicitly  write  the  matrix  elements  in  a simple  form. 
These  two  aids  to  the  organization  of  the  finite  element 
method  have  resulted  in  a saving  by  a factor  as  large  as 
twenty  in  the  computation  time  over  that  required  by  the 
method  described  in  [7],  even  though  the  two  approaches  can 
be  verified  to  be  analytically  and  numerically  identical. 


2, 4 A Comparison  of  Numerical  Methods  for  Application  to 

Miss  1 le  Plume  Problems 

In  this  section  we  take  a closer  look  at  some  of  the 
features  of  the  various  methods  indicated  earlier,  as  they 
apply  to  the  specific  application  we  have  in  mind,  namely, 
the  calculation  of  the  current  indued  on  a missile  with  an 
attached  electrically  inhomogeneous  exhaust  plume. 

We  may  immediately  rule  out  the  use  of  the  extended 
boundary  condition  or  T-matrix  approach  [2]  for  two  reasons. 
First,  the  method  is  slowly  convergent  when  the  scatterer 
is  not  nearly  spherical  and  hence  is  not  suitable  for  ap- 
plication to  the  thin  missile/plume  configuration.  Secondly, 
the  entire  domain  basis  representation  of  the  fields  (typ- 
ically in  terms  of  spherical  harmonics)  used  in  the  method 
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is  also  numerically  unstable  whenever  there  are  geometrical 
singularities  such  as  would  occur  near  the  tip  of  the  mis- 
sile and  at  the  missile-plume  junction  point  near  the  rocket 
nozzle. 

By  contrast  to  the  T-matrix  approach,  the  volume  polari- 
zation current  {-pproach  [1]  is  a very  numerically  stable  ap- 
proach. However,  it  generally  requires  a large  number  of 
unknowns  since  the  method  requires  one  to  solve  for  all 
three  components  of  an  effective  polarization  current.  In 
the  missile/plume  problem,  which  can  be  treated  as  a body 
of  revolution,  if  only  the  circumferentially  uniform  Fourier 
component  of  the  missile  current  is  desired,  the  number  of 
vector  components  in  the  polarization  current  is  reduced 
to  two.  Nevertheless,  several  other  factors  weigh  heavily 
against  this  approach.  The  first  is  that  the  matrix  that 
must  be  solved  is  full  (i.e.  not  sparse).  Not  only  does 
this  fact  means  that  matrix  fill  time  becomes  expensive,  but 
also  because  of  the  large  storage  requirement,  ou^-of-core 
matrix  solution  techniques  would  be  necessary.  Secondly, 
the  density  of  points  at  which  the  polarization  current  must 
be  sampled  is  related  to  the  local  wavelength  and  skin 
depth  in  the  medium  and  the  rate  at  which  the  local  medium 
parameters  are  changing.  Thus  the  method  is  not  suitable 
for  layered  inhomogeneous  bodies  or  regions  where  the 
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parameters  vary  slow  enough  that  layers  can  be  used  to 
approximate  the  scatterer. 

This  leaves  for  our  consideration  the  unincment  method, 
discussed  in  detail  in  the  previous  section,  and  the  surface 
integral  equation  formulation,  both  of  which  are  numerically 
stauie  and  result  in  banded  matrices.  Before  comparing  the 
methods,  we  consider  the  surface  formulation  in  further 
detail.  The  approach  for  layered  inhomogeneities  is  closely 
related  to  the  surface  field  approach  for  homogeneous  lossy 
dielectric  scatterers.  For  homogeneous  scatterers,  the 
approach  proceeds  as  follows.  Referring  to  Fig.  2.8,  one 
postulates  surface  currents  J and  K on  the  surface  S.  By 
relating  the  electric  and  magnetic  fields  to  J and  M and 
applying  the  boundary  condition  on  S that  the  tangential 
fields  have  to  be  continuous  across  S,  one  obtains  a pair  of 
coupled  Integral  equations  in  J and  M.  Using  numerical 
methods,  one  may  then  solve  the  coupled  integral  equations 
for  J and  M.  To  extend  the  approach  to  layered  inhomogeneous 
scatterers,  such  as  the  five-layer  scatterer  shown  in  Fig. 
2.9,  one  may  simultaneously  solve  for  postulated  surface 
currents  on  each  of  the  surfaces,  using  the  coupled  integral 
equations  derived  by  matching  the  fields  at  all  the  inter- 
faces. The  dimensions  of  the  resultant  matrix  are  such  that 
a full  matrix  of  those  dimensions  would  quickly  exceed  the 
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Figure  2.8.  A homogeneous  scatterer. 
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storage  available  on  most  computers  for  most  problems.  How- 
ever, if  one  recognizes  the  fact  that  a body  with  a layered 
inhomogeneity  has  a moment  matrix  which  is  block  tridiagonal, 
one  needs  only  to  simultaneously  store  three  of  the  much 
smaller  blocks  making  up  the  matrix  [4],  To  illustrate  this, 
we  consider  in  some  detail  the  example  of  the  five-layered 
inhomogeneous  scatterer  shown  in  Fig.  2.9.  If  is  the 

moment  matrix  for  sources  on  layer  j and  field  points  on 
layer  i,  then  the  overall  moment  matrix  equation  is  of  the 
form 
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(2.33) 


where  1^  is  the  column  vector  of  unknown  expansion  coef- 
ficients of  the  postulated  currents  Mi  on  surface  SA  and 

Vinc  is  related  to  the  tangential  incident  field  on  the 
outermost  layer.  We  note  that  the  matrix  is  block  tridi- 
gonal  in  form.  Beginning  with  the  equation  obtained  from 
the  first  row  of  (2.33),  repeated  elimination  results  in 
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*n  ^nn  ^n,n-l  ^n-1  ^n-l,n^  ^n,n+l  *n+l 

« Mn  Ir+1  , n * 1,2, 3, 4,  (2.34) 

with  M * 0,  Lc  , -I,  the  identity  matrix,  and  I . *-Vinc 

o 5,0  y n+1 

for  n*5.  We  note  immediately  tha'.  the  dimensions  of  the  matr- 
ices which  need  to  be  in  core  at  each  stage  of  the  iteration 
are  much  smaller  than  that  of  the  overall  matrix  in  (2.33). 

It  should  be  further  noted  that  the  iterative  approach  in 
(2.34)  uses  only  the  non-zero  sub-matrices  in  (2.33),  thus 
avoiding  storage  of  any  of  the  zero  sub-matrices  in  (2.33). 

In  order  to  compare  the  accuracy  and  efficiency  of  the 
surface  current  formulation  approach  with  the  univoment 
method,  a few  test  cases  were  tried.  Figs.  2.6  and  2.7  show 
the  scattered  fields  computed  by  the  two  methods  for  various 
two  dimensional  objects.  Based  on  our  experience  of  testing 
both  approaches,  we  offer  a list  of  observations  presented 
in  Table  2.1  as  a guide  to  choosing  between  the  two  methods. 

We  note  that  a number  of  variations  on  the  approaches 
considered  here  are  possible  and  these  can  significantly 
affect  our  conclusions.  For  example,  Morgan  [8]  uses  a 
homogeneous  core  region  in  his  unimoment  method  application 
and  obtains  a significant  saving  in  the  number  of  unknowns 
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required.  Some  of  the  features  of  the  various  methods  dis- 
cussed here  can  ai  . „ombined.  For  example,  one  might 
use  a finite  element  mesh. as  in  the  unimoment  method,  to 
determine  the  fields  interior  to  the  inhomogeneous  region 
and  use  a surface  current  formulation  to  simultaneously  de- 
termine the  fields  at  the  boundary  of  the  region  and  to 
enforce  the  radiation  condition.  Such  an  approach  seems  to 
combine  many  of  the  most  desirable  features  of  all  the  meth- 
ods discussed  here,  but, unfortunately , has  never  been 
tested*.' 

Of  the  two  methods  compared  in  Table  2.1,  we  choose  to 
use  the  surface  integral  equation  approach  to  treat  the 
missile/plume  problem.  This  choice  is  made  because  the 
missile/plume  configuration  is  a thin  structure  and  there- 
fore nor.  so  suitable  for  the  unimome.it  approach.  Further- 
more, we  are  primarily  interested  in  surface  currents  and 
not  with  scattered  fields.  For  these  reasons,  together  with 
some  of  the  complexities  associated  with  extending  the 
unimoment  approach  to  bodies  of  revolution,  we  have  chosen 
to  employ  the  surface  current  formulation  in  the  following 
chapters,  which  formulate  the  missile/plurae  problem. 
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CHAPTER  III 

NUMERICAL  SOLUTION  PROCEDURE  FOR  SCATTERING 
BY  LAYERED  DIELECTRIC  BODY  OF  REVOLUTION 

In  this  chapter,  the  integral  equations  are  formulated 
for  the  surface  currents  induced  by  a plane  wave  incident 
on  a dielectric  layered  body  of  revolution.  The  equivalence 
formulation  results  in  interactions  only  between  adjacent 
surfaces  and  hence  we  formulate  the  interactions  for  only 
two  such  surfaces,  the  pth  and  (p-l)th.  Many  of  the  details 
of  the  numerical  solution  procedure  essentially  parallel 
those  described  in  (llj  and  hence  are  not  explained  in  de- 
tail. Numerical  results  for  various  configurations  of 
layered  bodies  are  obtained  and  compared  with  other  results 
where  available. 

3 . 1 Formulation  of  the  Integral  Equations 

Fig.  3.1  shows  two  layers  of  a N-layered  dielectric 
body  of  revolution.  Thesu  layers  are  formed  by  revolving 
the  generator  arcs  ABC  and  PQR  about  the  z-axis  (assumed  to 
be  the  axis  of  the  body  of  revolution).  The  regions  bounded 
by  surfaces  S and  S , are  labelled  as  regions  p-1*  Pi  p+1 

p p-1 

going  from  the  inner  to  the  outer  layer.  Note  that  there 


p 

Figure  3.1.  Two  layers  of  a N - layered  dielectric 
body  of  revolution. 
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may  be  other  regions  interior  to  p-1  and/or  exterior  to  p+1. 
The  physical  parameters  of  the  regions  shown  are  e^,  o^, 

i = p-1,  p,  p+l . 

Boundary  conditions  require  that  the  total  fields  tan- 
gential to  the  two  interface  surfaces  be  continuous.  This 
implies  that 
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where  (E^,  H^),  i - p-1,  p,  p+1  are  the  total  fields  inside 

A A 

the  respective  regions,  and  n^  - <f>^  * t^  is  the  outward  unit 

t h a 

vector  normal  to  the  i surface,  with  t^  being  the  unit 
vector  along  the  generator  arc.  On  the  N-th  layer,  we  have 
similarly, 


ftN  X (1N  " EN+1)  * fiN  * * 


inc 
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Using  the  equivalence  principle,  we  divide  the  problem  into 
three  parts,  one  each  for  the  evaluation  of  the  fields  in 
the  three  regions.  Thus,  if  we  assume  equivalent  electric 
and  magnetic  currents  , i = p-1,  p,  on  the  two  sur- 

faces, the  three  equivalences  are  obtained  as  shown  in  Fig. 
3.2.  Note  that  the  currents  all  radiate  in  a homogeneous 
medium  in  each  of  the  cases  depicted  in  Fig.  3.2.  The  equi- 
valent currents  on  the  interior  of  a given  surface  are  merely 
the  negative  of  those  on  the  exterior  to  the  surface  in  view 
of  the  continuity  of  the  tangential  fields  at  the  surfaces. 
The  fields  in  (3.1)  and  (3.2)  can  thus  be  written  as 


Et(?)  - 

-juAi(r) 

- ?♦*(?) 

1 

" ei 

VxFi(r) 

(3.4a) 

Hi(r)  - 

-jo)Fi(r) 

- V4»“(r) 

+ -±- 

VxAi (r) 

(3.4b) 

i - p-1,  p,  p+1,  (3.4b) 


where  the  potentials  are  defined  as  follows: 


V1. 


S 

P 


1 


G .(i,  *')  ds' 


(3.5a) 


.-r«w'w,vi 
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// 


FP-i(')  ’ - ~fr  J J > d8' 

Vi 


$ (r)  = - 

p-lv  1 4tre 


' If  ' 


pp-lfr,)Gp-l(?»;  > ds’ 


P-1 


Cl(_r) 


4iry 


/; 


P p 

SP-1 


Pp-1<;  )GP-1(?,?  > d8’ 


V*}  *"//  Vi(i,)V**?,)  d8‘ 

;p-l 


// 


- ^M3p<f,)Gp(5-f,)  ds' 


% (r)'4?jj  flp-l('')Cp<f'f'>ds' 

Vi 


// 


I I ds’ 


(3.5b) 


(3.5c) 


( 3 . 5d) 


(3.6a) 


(3.6b) 
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f 


pp-l^  ' ^ds  ’ 

Vi 


Pi 


4 *kJJ  p^<?,)  gp<?,?,)  ds 

s 

p 


*lCr)  JJ  pp-i(?,)fiP(?>?,>  d*’ 


//• 


1p“(r')G„(r,r,)ds' 

P.  ' 

S 

P 


47fP_  IIP  ' ~p 


// 


’p+l"'  4tt  II  Jp(r,)Gp+i^.?')d8' 


f. 


Vi<;>  ■ ~4¥ii/Rp(;’>Gp+1<;.I')<i»' 

S 


♦“,<;>  - 1 * * e-  - 


//■ 


P+1"'  47rep  + I / fPp(r’r,)Gp+l(?’^)ds 


<■",(?>  = 1 


//. 


P+l'w  " 4 7ipp+l  JJ  Pp(r’r  ’ )Gu+1(r»r ’)ds 


(3.6c) 


( 3 . 6d) 


(3.7a) 


(3.7b) 


(3.7c) 


( 3 . 7d) 
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and  where 


-jt^R 

Gi(?,?’)  = (3.8) 

i * P-l»  P»  P+1 

Sj 

R = |r-r'|  = [p2  + p'2  - 2pp*  cos  + (z-z*)2]  . 

Using  the  continuity  equation,  one  can  express  the  charge 
densities  in  terms  of  the  current  densities  as 

Pi  " i ^'a  * V*’)}  (3.9a) 

p"  - - [V  • Mi(r’)]  (3.9b) 

1 w S * 

i * P-1,  P. 


Using  (3.4)  through  (3.9)  in  (3,1)  and  (3,2),  one  obtains 
the  integro-dif ferential  equations  for  the  unknown  electric 
and  magnetic  currents,  shown  on  the  following  pages. 


<e 
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j 1 to 


If 


Vi  ' \iiJJ  jh(Vi  Vi  + Vp)ds 

sP-i 


+ 4^3  VI  I (?;  • Vl> 


•// 

Vi 

■// 


G i G , 

+ d3' 
P-1  P, 


fip«=p.i + S)ds’  - tijJWV’’ 


IS 


P-1 


4ttui 


// 


V I I (VI  * J„)  7^  ds  ' “ 7 — V * I I M G ds 

8 P ep  4”  JJ  P P 


If 


0,  reS 


P-1 


(3.10a) 


>-l  * //  Sp-l(ep-l  Gp-1  + 


e G )ds' 
p p 


P-1 


(3.10b) 
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n x 
P 


Ijw 

I 4u 


I I J (M  G + (J 

JJ  ? V V 


p+i  GP+i)  ds' 


G G . - 
_£  + _E±i  ds. 

P P+1j 


+ — V I I < V * • J 

4ttu)  J J s P 
s 

P 

+ t7  7 x//fiP(GP  + °p+i)ds'  ■ iVp  ds' 


P-1 


4ttu) 


//<*;  • vi>  i2-*8'  - x JJ  Vi 


G ds 
P 


*P-1 


P-1 


reS  (3 
P 


f. 


II 


n * I I M (c  G + e . - G .-Jds' 

P |4tt  JJ  pv  p p P+1  P+1 
1 S 

P 


+ - J-  V I I (V*  • M 

4ttu)  «/ ,/  s I 

S 

P 

//> 


G G . 

+ .£.tl 

^p  ^ p + 1 


ds  ' 


7-  V x | |J  (G  + G ,.)ds’  - 
4TT  JJ  p p P+1 


ico 


// 


// 


, . . M . c G ds' 

4tt  JJ  p-1  p p 

!P-1 


(V  • M ,)  ds 
s p-1  W 


4nw  JJ  ' s p-1'  U 
‘p-1 


* r»  v * //  Vi 


P-1 


G ds 
P 


reS  . ( 3 . 
P 


11a) 


■|=0, 

11b) 
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In  (3.10)  and  (3.11),  the  dependence  of  the  various  quanti- 


ties on  source  and/or  field  coordinates  is  implicitly  under- 
stood. We  next  note  that,  for  i ^ j, 


//■ 


V xj,  UjG(r.,rj)ds* 


// 


Uj  x \,G(ri,rj)ds' 


(3.12a) 


Uj  = Jj  °r  MJ* 


wherein  the  curl  has  been  taken  inside  the  integral  since 
the  field  point  r^  and  source  point  r^  are  on  different 
surfaces  and  hence  G(r^,rj)  is  non-singular,  assuming  that 
the  layers  do  not  touch  one  another.  However  when  r and  r' 
are  on  the  same  surface,  then  it  can  be  shown  [12]  that  on 
the  surface  S, 


//[ 


5 * JJ  Vci+<Wds' 


Ui  x V(Gi+G1  + 1)ds’  (3.12b) 


where 


indicates  a deleted  integral  around  r.  Equa- 


tions (3.10)  and  (3.11)  can  oe  written  in  compact  form  as 


S , .( 3 , , M 1)+C.(J,M)“0  (3.13a) 

p-l.p-l  P-1  P-1  P-l.P  P p 


'p,p-l(j  .,  M .)  + S (J  , M ) = 0 (3.13b) 

p-1  p-1  p,p  p p 


i 
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where  S ^ ^ ( J ^ , M^)  represents  the  contribution  to  both  elec- 
tric and  magnetic  fields  on  surface  i from  currents  existing 
on  surface  i,  while  0^(3^,  M ^ ) represents  the  electric  and 
magnetic  fields  evaluated  on  surface  i due  to  currents 
existing  on  surface  j.  Such  a compact  form  is  necessary  in 
order  to  generalize  the  procedure  to  an  arbitrary  number  of 
layers.  We  can  write  (3.10)  and  (3.11)  in  component  operator 
form  as  follows: 


-V 


cij<3j 


s 


s 


s 


s 


n(Jit> 

s12(Jt> 

<J> 

S22(J;j  ) 

31  <"*  i t > 

S32(Ji  * 
0 

*.<V 

S42(Ji  > 
$ 

n(V 

Cl2(Jj  5 

21<Jj,) 

C22(J , ) 
♦ 

31UJt) 

C J2(Jj J 

41(V 

c42  ^ J j ^ 

<P 

13 (Mi t ^ 

s14<H1  > 

<P 

23(Mi^) 

s24(M1  > 

<P 

33(Mit) 

s 3 4 ( M i J 
<P 

43(Mij.) 

S44(Mi  > 

<t>  . 

(3.14) 

13<V 

C“(V 

2 3 (M j * 

J t 

C»°V 

33(Hjt) 

C»‘V 

*3°V 

‘"‘V 

(3.15) 
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with  i = j+1  or  j-1.  Herein  the  first  and  third  row  of 
operators  yield  the  t-directed  component  of»the  tangential 
electric  and  magnetic  fields,  respectively,  while  the 
second  and  fourth  row  of  operators  yield  the  ^-component  of 
the  tangential  electric  and  magnetic  fields.  The  dependence 
of  the  operators  and  on  the  various  current  compon- 

ents is  clearly  indicated.  We  next  express  all  the  operators 
in  (3.14)  and  (3.15)  in  terms  of  the  local  coordinates  ( t , ) 
on  the  surfaces.  An  orthogonal  system  of  unit  vectors 
(n^  t^),  i * p-1,  P,  are  associated  with  each  point 

(t^,  * * p-1,  p on  the  surface  S^.  These  unit  vectors 

are  defined  as  follows: 


n^  = cos  y^cos  $ 5c  + cos  y^sin  4>i  y - sin  y^  £ (3.16a) 

= -sin  4>j  x + cos  y (3.16b) 

* sin  Yicos  & + sin  YjSin  4*^  y + cos  £ (3.16c) 

i = p-1,  p,  where  y^  if*  the  angle  between  the  tangent  to  the 
generating  arc,  and  the  z-axis,  defined  to  be  oositive 

a A 

if  t ^ points  away  from  the  z-axis  and  negative  if  t ^ points 
toward  the  z-axis.  The  surface  divergence  can  now  be  writ- 


ten as 
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v • n 

si  1 


i.  _9_  . , „ v , 1.  9 

pj  at;  (»i  V p;  a*, 


r <Ux  ) 


(3.17) 


U1  = J or  Mi,  i 


P-1.P 


Exp-mding  (3.10)  and  (3.11)  into  component  form  and  compar- 
ing with  (3.14)  and  (3.15),  one  can  easily  obtain  the 
expressions  for  s.  and  c . For  notational  simplicity,  we 
introduce  the  following  operators: 


g11(U;  p , e , tif  t’,  4^,  <)»p 


II 


kaJJV  lSin  7i8ln  Yj  cos(4»  ) 


+ cos  71cos  Yj ] pG  ( t ^ , t j ) ds ' 


4ii(d  3t 


II 


i a 


r (p;  U) 


PJ  9tj  j 


ds’ 


(3.18a) 
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8 , 2 (U ; y,  e,  t*  , ^,4*]) 

‘ - 4JS//U  sin  Y1sin(4>j— 4)i>  liG(t1(  t * )ds  1 


4tt 


’j 


// 


+ Tnw  3tt  JJ  P]  34>]  (U)  e 

Sj 


G(t^-  ds’ 


8n(U;  h,  e,  ti,  i‘,  4^.  4>]) 


WJ  ft(pj  “ln  YiC°“  Yi  ' P‘C°8  YiSl"  Yj> 

s 


j 


x + (z1-zpsin  Y^in  Y j s in  (4>'j-'P  £)  ] 


dR 


14 


(U;  ;j  , c > t ^ , t , ^ j ) 


4tfJT  R 


~ (pj  cos  ya  " P±cos  Y^os^j  -<i) 


f (zi-*z’)sin  Y i cos(<})j  - <J>  ± ) ] ^ G(ti,tj)dsl 


■ . w i irr -m~~**xr*** 


(3.18b) 


(3.18c) 


(3. 18d) 
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821(U;  ii,  e,  tj,  4>± * <J>j) 


= 1“ 


//■ 


SJJ  U sin  y! sin^j-if^)  pG(ti>tj)  ds' 
S 


j 


+ d^ffp^  3tj  (pju) 

Sj 


G(tt,tp 


ds' 


B22(U;  u,  e,  tif  t',  ♦i.  *') 


// 


i^JJv  coat*"-*,)  UG<t1,tpds’ 

SJ 


+ i 


Antop^  3 <|>  ^ 


f f — — (u) 

' 

JJ  oj  »♦! 

e 

, 4 

ds  ' 


023(U;  m , e,  ti,  tj,  4>i,  $’) 


if- 


«„ jj  r i-v”  - "J'08  \;'=‘>8<*;-*1) 

s 


j 


- (zt-zp  sin  YjCos(4>j-4>1)  ] ^ G(ti»tj)ds 


( 3 . 18e ) 


(3 . 18f ) 


I 


( 3 . 18g) 
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fi24(U;  y,  c,  tjL,  t’,  <j>i,  <j>p 


Gtt^tpds'  ( 3, 18h) 

’3 


where 


-jkR 


G(ti*ti)  " R 


(3.19a) 

1 


R = lri"rj 


[p2  + Pj2  - 2PtPj  cos  (♦1-4> ! ) + (zi-  *j)2]2. 

(3.19b) 


With  the  above  definitions  we  have  the  various  and 

as : 


sn(JP-it)  " 6u(JP-iti  MP-i’  eP-i’  Vr  t'P-i’  Vi’  ^ ' p-i) 

+ B,.(J  , ; u , e , t , t'  1f  4>  ♦'  .) 

ll  p-it  P P P_  1 p- 1 p-i  p-i 


(3.20a) 


12 


(J  I ) * Bj2^r_i  ’ M_  i > ^n-1’  ^o—l*  ^ o-l*  ^o-l*  ^ n-1^ 

p-  14  p—1  P”-1  P“  1 r — 1 P~  1 r~ i ‘ 


f 6.  »(l  . » M > £ » t t « » 4*  i I 4*  l) 

12  p-1^  p p p-1  p-1  p-1  p-1 

(3.20b) 
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!13(Mp-lt)  " ei3(MP-lt;  Pp-1’  £p-l’  tp-l’  t'p-l’  *p-l*  ‘Vl) 

+ 0 . - (M  . y , e , t t ' , <j>  4>'  , ) 

13  P_1ti  P P P-1  P"1  P~ 1 P 1 

(3.20c) 


■ ‘u’vv  vi’  Vi-  ‘p-i-  ''p-i-  Vi-  ♦p-i> 


14  p 


<*» 


+ Bu(Mp_1  ; yp,  sp,  tp  ^ t'p_1,  p_ x , 4»'p_1) 


(3 . 20d) 


921(Jp-l  * " 621(Jp-l  5 Up-1’  ep-l’  tp~l,  t'p-l>  %-!’  ♦'p-l* 


+ 601(J„  , . y . e . t ,,  t ' , 6 , <j>'  .) 

21  P-1t>  P P P-x  P~ 1 p-1  P-1 

(3. 20e) 

s.,_(J  . ) - , *,  y„  e , t , t'  , 4>  <J> ' .) 

22  p-1  22  p-l$  p-1  p-1  p-1  p-1  p-1  p-1 


+ P22(Jp-l^;  V %’  %-!’  t'p-l’  ^p-1’  <J,'p-l) 

( 3 . 2 Of  ) 


S23(Mp-lt)  = 823<Mp-lt:  up-l’  Cp-1’  tp-l’  ^p-l’  ^p-1*  '?>,p-l) 

+ *23<Vv  V V Vi*  "Vi*  Vi- 


( 3 . 2 Og) 


62 


S43(Mp-lt)  = e21(Mp-lt;  £p-l»  V'p-1*  tp-l*  t'p-l*  ♦p-l*  <{>’p-l) 


+ 621(MP-lt;  V V tp-l*  t,p-l*  Vl’  ^’p-^ 

(3 . 20o) 


S44(Mp-1(J))  “ 622(Mp-l(j);  £p-l»  Up-l*  tp-l’  ^p-l’  *p-l’  <t>’p-l) 


+ 622(Mp-l(J):  £p*  V tp-l»  t'p-l»  ♦p-l*  ♦’p-l) 


( 3 . 20p) 


CH(Jpt) 


"8u(Jpt;  v v Vl’  cp’  Vi*  (3-21a) 


C»(V 


6.»(J  > P > E > t . , t ' , <J>  . i <$>  _ ) (3.21b) 

12  P<t>  P P P_1  P P"1  P 


cn(M  ) 

1 J pt 


6i3<Hpt!  V V Vi’  t’,‘  Vi-  *p)  (3-21c) 


c..(M  ) - -0  . (M  ; M , e , t t',  4>  <J>')  , - „ . , s 

14  p^  14  p^  p’  p p-1  p’  Yp-1  rp  ( 3 . 2 Id) 


c-.(J  ) 

21  Pt 


-8  (J  ; y , e , t t\  <b  , <J>’)  (3.21e) 

21  Pt  p p p-1  p p-1  P 


22  p^  2 ” P (p  P P P"1  P P'1  P 


C23(Mp  ) 
Kt 


623(Mpt;  %’  ep’  ‘p-l*  tp*  %-!’  (3,21g) 
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cn(H  ) = 

24  p, 

-624<Mp^  %• 

V 

Vl' 

V 

Vi-  V 

( 3 . 21h) 

c31<Jpt>  ■ 

’ C13(Jpt) 

(3.211) 

C.^J  ) = 

41  p* 

- C, . (J  ) 

14  p* 

( 3 . 21  j ) 

33("pt) 

■6n(Mpt;  V 

P’ 

'p-1' 

V 

Vi-  v 

(3.21k) 

24<V  " 

-'»(V  v 

V 

Vi' 

v 

Vi’  V 

(3.21il) 

-C23(Jpt) 

(3.21m) 

V<Jp  > ‘ 

<P 

( 3 . 21n) 

_e21(“pt*  V 

V 

Vi' 

v 

Vi'  V 

(3. 21o) 

:44<V  ■ 

‘V<Mp/  V 

V 

Vi' 

*;• 

Vi'  ♦;>• 

(3 . 21p) 

We  note  (3.20)  and  (3.21)  define  S . ( J .,  M ) and 

p—l^p—l  p — 1 p — 1 

C . , (J  , M ) in  (3.13a).  To  find  similar  expressions  for 

p-1  p p p 

S (J  , M ) in  (3.13b),  we  replace  (p-1)  by  p and  p by  (p+1) 
P » P P P 

in  (3.19).  However,  we  replace  (p-1)  by  p and  p by  (p-1)  in 


(3.20)  only  for  the  field  and  source  coordinate  variables, 
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maintaining  the  medium  parameters  to  be  e^,  to  obtain 

C , (J  . , M .)  in  (3.13b) . 

p,p-l  p-1  p-1 

We  next  turn  our  attention  to  the  numerical  evaluation 
of  the  various  quantities  in  (3.14)  and  (3.15).  For  this, 
we  divide  the  generating  arcs  into  a number  of  linear  seg- 
ments as  shown  in  Fig.  3.3.  The  points  t^  , t.  , . ...t^ 

o 1 n 

...,  i ■ p-l,p  specify  the  end  points  of  the  linear  segments. 
The  half-points,  t.  , t t,  ....,  i = p-l,p,  are 

\ n+h 

defined  by 


2 


i 

i 


p-1,  1 < n < N .+1 

" - P"1  (3.22). 

p,  1 < n < N +1 
- P 


We  next  define  expansion  functions  for  the  electric  and  mag- 
netic currents.  Since  the  scatterer  is  a body  of  revolution, 
we  choose  the  $-  variation  of  the  currents  to  be  the  Fourier  modes 
in  <|> , while  the  t-variation  is  expanded  in  a pulse 
basis  set.  Thus  we  have 


V‘i-  *p  ■ 2^  2 2 


m “ n ■ 1 


1™  pn(  ejm^ 
lt  1 1 


+ $ 


„ N.  + l 

. I 2 I 

n»_oo  n*l 


j"n  p"(t;>  eJ"+; 


(3.23) 
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Figure  3.3.  Linear  segmentation  of  the  generating  arcs. 
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where  N,  = N - or  N for  i = p-1,  p.  The  derivative  of 
i P-1  P 

J.  with  respect  to  t is  approximated  as 
1t  1 


N £+l 


afx'o;  %<';•  +;>)  - h 2 2 


m = -°° 


n - 1 


..inn  .m,n-l 

i 1 i 

t t 


At 


(3.24) 


with  I®*°  « j“*Ni+l  * 0.  Herein  P^(t!)  and  P^(t')  are 
it  it  1 i 2 i 

defined  as  follows,  for  i * p~l,p: 


FI<*i> 


i .i'ii'i 

n-»i  n+v 


1.  t 

0,  otherwise 


(3. 25a) 


P2(tI> 


. < t’  < t 

i — i - n 

n-1 


1,  t 

0,  otherwise 


(3.25b) 


with  At. 

l 


Ci  - Ci  ' 
n n-1 


■ Hp,  - p,  )2  + («,  - «,  >2l\ 

n n-1  n n-1 


(3.26) 


lmn  and  Jp  are  the  unknown  expansion  coefficients  of  the 
1 1 % 

"total  modal  current"  (2it  J ) in  the  t-direction  and 

n 1 1 
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Che  actual  modal  current  density  in  the  ^-direction.  The 

expansion  of  t-  and  ^-components  in  terms  of  the  "total" 

current  and  current  density,  respectively,  arises  as  a natural 

choice,  particularly  when  one  applies  the  surface  divergence 

to  the  current  representation.  One  observes  further  that 

J,  (t!,  <f>!)  , i * tj>-l,p,  appears  in  all  integrals  with  an 
Xt  1 1 

associated  factor  p^,  thus  making  it  natural  to  consider  the 

product  p^J^  (tj,  4*1)  the  unknown  current.  Ji  (t^,  $[)  > 
t 4> 

on  the  otherhand,  does  not  appear  this  way  and  the  artificial 

introduction  of  factor  would  introduce  unwanted  singulari- 
ties in  some  of  the  integrals.  Hence  we  expand  J alone. 

The  magnetic  currents  M^(t^,  and  their  derivatives  are 

expanded  similarly  with  K°n  * 2ir p * M°n  and  m"111  replacing 

t 1 t i<(> 

l”n  and  J®n  respectively. 

t S 

Features  of  the  current  expansion  scheme  in  (3.23)  and 
(j,24)  are  further  discussed  in  [11].  The  salient  points 
are  that  since  the  Fourier  modes  are  decoupled,  one  can 
solve  for  the  Fourier  components  mode  by  mode.  The  staggered 
pulse  basis  set  in  the  t-direction  permits  the  modeling  of 
open  and  closed  (conducting)  bodies  and  bodies  with  sharp 
edges  without  the  placement  of  observation  points  for  poten- 
tial quantities  at  points  wherein  the  corresponding  sources 
(current  or  charge)  may  be  singular. 
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We  next  define  the  testing  functions 

- j p<J> 

T^q  = Pq(t  )e  j (3.27a) 

J ^ ■**  J 

— i p d) 

T?q  - Pq(t  )e  i . (3.27b) 

J2  1 J 

The  t-c.  ooponent  equations  in  (3.13a)  and  (3.13b)  are  tested 
with  (3.27a)  while  the  ^-component  equations  in  (3.13a)  and 
(3.13b)  are  tested  with  (3.27b).  In  order  to  perform  the 
integrations  analytically,  we  expand  the  Green's  functions 
in  Fourier  series,  as 
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(3.28a) 
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- 1 k . R 
1 o 


m ' 


/I  d [ i.  o\ 

R~"  dT'  |~R 1 cos  (3.28b) 

o o'  o ' 


Ro  = l Pi  + pj  - 2pipj  cos  ^zi“zp2]^  . 


(3.28c) 


where  £ = <}>^-  4> ^ . The  form  of  the  integrand  in  the  above  is 
obtained  by  noting  that  exp (- j kiRQ) /Rq  is  even  in  £.  With 
this,  the  surface  integrals  now  reduce  to  a single  integral 
over  t.  Following  the  procedures  similar  to  those  used  for 
a homogeneous  body  of  revolution  [11],  we  can  write  the 
expressions  for  the  elements  flj?  , where  m refers  to  the  m1"*1 
Fourier  mode,  while  q and  n refer  to  the  field  and  source 
indices  respectively.  For  notational  simplicity,  we  define 

f j2 

«/•  i < t . * t ; t ; a)  -I  Gi  (t  , t * ) d t * (3.29a) 

J1  J2  o m q J J 


<P,  (t,  , t ; t ; m) 
Ji  J2  q 


r 2 

J Gt  (tt  , dtj  (3.29b) 


m q 


where  (t  , t')  are  defined  in  (3.28).  We  introduce  next 
m 3 

the  weight  functions  arising  due  to  the  testing  in  the  t- 
direction. 
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Xs<Ati  ’ ) 

q q 


At  ^ sin  + Atj  sin  ^ 

a +i _3 g a 


(3.30a) 


XcUti  • Yi  > 
q q 


At.  cos  v.  + At.  sin  y. 
i , , ' i i ' i 

q*1 d 9 


a. . 


(3.30b) 


The  various  can  now  be  written  explicitly.  One  notes 

that  the  0^*j  defined  in  the  following  pages  are  of  the 

IQ 

same  form  as  the  0**”  defined  in  (5.25)  of  [ 1 1 ] , except  that 

■*  m 

the  elements  there  contain  contributions  from  both  sides  of 
a surface  and  therefore  involve  two  medium  properties.  The 
reason  for  defining  the  0^”  in  terms  of  the  parameters  of 

J IQ 

a single  medium  is  that  s ^ involves  the  computation  of  0 
for  two  media  , while  c.^  involves  the  computation  of  0^ 
for  only  one  medium.  Thus  the  same  computer  subroutine  for 


m 


computing  0^  can  be  used  either  for  the  computation  of 
s or  c...  Hence  we  finally  write  the  0^"  as  follows: 
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$2"  = 1 ± [y^iJ(t  ,t  ; t£  ; m+1) 
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024  * " I jAti  Ul1<ti  » 5 S m>»  (3 . 31h) 

m q Jn-1  Jn  q— *s 

q * 1,2.. . N^+l 
n » 1 , 2 . . . N j+1  . 


The  U's  in  (3.31)  are  defined  aa  in  (5.26)  of  [11]  with 

(t  , t ) replaced  by  (t  , t ),  t replaced  by  t.  , the 
1 * •5l;52q  q 

field  points  (p,$,z)  replaced  by  (Pi»^i*zi)»  the  source 

points  (p » , z * ) replaced  by  (pj,  <J>j  , z ! ) and  (5.27)  of 

[11]  i 8 replaced  by 


G 


(3.32) 


where  k corresponds  to  the  medium  p-1,  p,  or  p+1, as  the  case 
may  be.  This  completes  the  definition  of  the  elements  of 
the  self- Impedance  and  mutual  Impedance  matrices. 

The  incident  plane  wave  fields  are  incident  on  the 


or  the  outermost  Interface.  These  can  be  expressed  as, 

(3.33a) 
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. i £i  i ?i . ”jkN+ln  * rN 
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HinC 
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(3.33b) 


where 
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0*  3 cos  8*  cos  (j>*  8 + cos  0*  sin  4>*  y - sin  0*  z (3.34a) 

<f>  = -sin  4>  x + cos  4>  y (3.34b) 

n = -sin  0 * cos  <J>*  x - sin  0*  sin  <f>*  y - cos  O'*-  z (3.34c) 

*rN  " PN  COS  *W  X + PN  Sln  y + ZN  Z ( 3 . 34d ) 

and  n„.,  is  the  characteristic  wave  impedance  of  the  (N+l) 
or  outermost  region.  The  tangential  components  of  (3.33) 
are  the  forcing  functions  of(3.,3).  The  corresponding  elements 
of  the  drive  vector  are  given  by  (5.30)  of  [ll],  where  the 
field  coordinates  refer  to  the  or  outermost  interface. 

After  computing  the  matrix  elements,  the  iteration 
procedure  given  in  (2.34)  can  now  be  used  to  solve  for  the 
currents  on  the  or  outermost  surface.  If  one  requires 

the  currents  on  all  the  1.. vers,  then  one  has  to  carry  out 
a back  substitution  process  beginning  with  the  currents  com- 
puted on  the  outermost  layer.  However,  if  one  is  interested 
only  in  the  fields  on  the  exterior  of  the  scatterer,  then  the 
currents  on  the  outermost  layer  are  sufficient. 

The  above  procedure  must  in  principle  be  carried  out 
for  each  of  the  modes  m,  £ m _<  °°.  In  practice,  however, 
one  achieves  convergence  with  only  a finite  number  of  modes. 
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Further,  as  shown  in  [ll],  computation  of  only  the  positive 
inodes  is  sufficient  to  obtain  the  total  contribution  from 
all  modes,  positive  and  negative.  This  is  because  of  vari- 
ous symmetry  relationships  between  the  matrix  elements  and 
currents  for  the  negative  modes  and  those  for  positive  modes. 


3 . 2 Eigenfunction  Solution  for  a Three-layered  Dielectric 

Sphere 

In  this  section,  we  obtain  expressions  for  the  induced 
currents  on  a three-layered  dielectric  sphere  illuminated  by 
a plane  wave.  This  is  done  to  obtain  a check  on  the  formula- 
tion and  computer  code  for  the  layered  body  of  revolution 
discussed  in  the  previous  section.  The  procedure  is  similar 
to  that  for  the  case  of  a homogeneous  sphere  [ 13]  which 
results  in  the  well-known  Mie  series  solution. 

Consider  the  inhomogeneous  sphere  shown  in  Fig.  3.4. 

The  various  media  are  characterized  by  parameters  (y^,  e^, 

0 ),  i = 1,4.  An  incident  plane  wave,  assumed  to  be  polari- 
zed in  the  x-direction  and  travelling  in  the  positive  z- 
direction,  can  be  expressed  as 
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where  is  the  characteristic  wave  impedance  of  medium  4. 
For  convenience  in  applying  boundary  conditions,  we  express 
the  fields  in  terms  of  TM  to  r and  TE  to  r problems.  Thus 
we  have  associated  incident  electric  and  magnetic  vector 
potentials,  F^,  A*,  respectively,  associated  with  the  inci- 
dent fields  in  (3.35)  as  [13], 
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(3.36a) 
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(3.36c) 


Here 


3n(kr) 


kr  j n (kr ) 


(3.37) 


where  j (x)  is  the  spherical  Bessel  function  of  order  n, 
n 

Jn(x)  is  the  ordinary  cylindrical  Bessel  function,  P^(cos  0) 
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is  the  associated  Legendre  polynomial  of  degree  n,  order 
one.  We  next  expand  the  fields  in  the  various  regions  in 
terms  of  potentials  of  the  same  form  as  those  expressing  the 
incident  field.  Thus  the  fields  in  the  various  regions  are 
expressed  in  terms  of  the  following  potentials: 
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Herein,  N^fkr)  and  (kr)  denote  kr  times  the  ordinary 

spherical  Neumann  and  spherical  Hankel  function  of  the 
second  kind,  respectively.  Note  that  the  fields  in  region  1 
are  rinite  while  the  scattered  fields  are  represented  in 
terms  of  outgoing  waves.  The  electric  and  magnetic  fields 
can  be  determined  from  the  potentials  as  follows  [ 13]  : 
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( 3 . 39d) 


The  bouudar'  conditions  to  be  met  are 
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O<0<TI,  — TT  < 4>  < TT , i = 1,2,3. 
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In  view  of  the  independence  of  the  TEr  and  TM^  modes,  we 
simply  require 
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i - 1,2,3  . 


Application  of  (3.41)  automatically  ensures  the  continuity 
of  the  (^-components  of  the  fields  as  well.  Use  of  (3.38) 
and  (3.39)  in  (3.41)  thus  yields 
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wherein  the  orthogonality  amongst  the  modes  has  been  utilized. 
Equafions  (3.42)  are  12  equations  in  the  12  unknown  coef- 
ficients. The  coefficients  were  numerically  determined  by 
solving  the  system  of  linear  equations  (3.42).  Jne  notes 
that  if  (3.42)  is  written  in  matrix  form,  the  matrix  can  be 
partitioned  into  two  independent  matrices  (involving  TEr  and 
TM^  coefficients)  which  are  individually  block  tridiagonal. 
Hence  an  iteration  scheme,  similar  to  the  one  used  for  the 
layered  body  of  revolution,  can  be  developed.  Equation 
(3.42)  is  to  be  solved  for  each  of  the  modes  n.  (Note  that 
only  the  0- dependence  cf  the  fields  varies  from  mode  to 
mode,  whereas  the  $ -dependence  of  the  fields  is  fixed. 
Because  we  have  considered  axial  incidence,  only  the 
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Fourier  mode  in  <t>  exists  as  indicated  in  (3.38)). 

The  number  of  modes  is  truncated  at  a finite  number  after 
convergence  has  been  achieved. 

Having  determined  the  coefficients,  we  can  compute  the 
surface  current  on  the  outermost  layer.  We  have 
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Using  (3.38)  and  (3.39)  in  (3.43),  we  get  the  various  com- 
ponents of  the  currents: 


E cos  6 

j m _j|  I m ® 

°3  ♦4  r“a3 


5!  (\3n(lV3> 


n * 1 


(2) 


+ rn  Hn  (k4a3}>  ^ P-(cos  0) 


E cos  <b 
o 

j top 3 sin  6 


d „1 


d 0 n 


I 

n *=  1 


a J ' (k  a - ) 
n n 4 3 


(3.44a) 


+ s 


H(2)l (k.a ,)}  P1(cos  6) 
n n 4 3 J n 


87 


r4s  * He 


r * a. 


+ r H(2)(k,a  ) 
n n 4 3' 


E sin  6 
o T 

sin  0 


P (cos  0) 
n 


00 

lu  i 

I n n 
n * 1 v 


(k4*3> 


E sin  6 
o \ 

j u>MAa3 


2.  (•■ 


<V3>  + % 


n - 1 


M - E. 
°3  *4 


r * a. 


E sin  <J> 
o 

k a 
13 


l , (•. 


(k^a3> 


n - 1 


+ 8n  9) 


E sin  4> 
o 

jkAa3  sin  0 


2 k 


( k 4 a ;. ) 


n * 1 


+ r fk2)  • (k 
n n 


4a3>}  PJ 


(cos  0) 


} ~h  pi(cos6) 

( 3 .'44b  ) 


(3.44c) 


88 


M.  * "EA 

‘t’A  9, 


E cos  d> 

0 T 

r.a-  sin  0 
4 3 


2 K s. 


(k4a3> 


s II (k  a ) 
sn  11  n U4a3; 


P (cos 
n 


6) 


E cos  I 
o _ 

jk4a3 


2 


n - 1 


a J'  (k.a  ) 
n n 4 3 


+ c„  ^2>,<V3>}  v>  pi<"»  e> 


( 3 . 44d) 


Equations  (3.44)  thus  give  the  surface  currents  on  the 
outermost  layer.  lhe  infinite  summation  is  truncated  to  a 
finite  number,  after  ascertaining  the  numerical  convergence. 

3 . 3 Numerical  Results 

The  numerical  procedure  described  in  Section  3.1  for  a 
layered  dielectric  body  of  revolution  has  been  incorporated 
into  a computer  code.  The  resultant  computer  code  is 
hereafter  referred  to  as  "LDBR".  We  compare  the  numerical 
results  from  the  code  LDBR  with  the  solution  for  the  layered 
sphere  developed  in  Section  3.2.  After  establishing  the  code, 
additional  geometries  for  whicn  exact  solutions  are  not 


available  are  considered. 
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The  various  scattering  models  considered  are  assumed  to 

be  illuminated  by  a uniform  plane  wave  travelling  in  the 

positive  z-direction.  Since  a plane  wave  incident  along  the 

axis  of  a body  of  revolution  excites  only  the  n=+l  Fourier 

+j$ 

mode  with  e variation,  the  currents  on  the  surface  have 
the  form  [11] 

J(t,<}>)  ■ J ( t ) cos  (#>■+•  j J^Ct)  sin  <f> 

M ( t , $)  * j Mt(t)  sin  <{>  + M^(t)  cos  4> 


where 

i 

J (t)  - 2 J®(t) 

P P 

, p - t or  , 

M ( t ) - 2 M®(t) 

P P j 

I 

0 0 

and  where  J and  M are  the  electric  and  magnetic  currents 
P P 

resulting  from  a 0-polarized  incident  field  with  a circum- 
ferential variation  of  the  form  e^.  The  figures  in  this 
section  show  Jt  and  in  the  $ = 0°  plane,  while  and 
M ar<‘  m the  4>  * 90°  plane. 

As  a check  on  the  accuracy  of  the  iterative  computa- 
tions, we  have  considered  some  cases  wherein  the  scatterer 
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is  made  up  of  dummy  layers.  That  is,  the  entire  scatterer, 
though  homogeneous,  is  treated  numerically  as  being  made  up 
of  layers,  each  layer  with  the  same  material  parameters. 
Fig.  3.5  indicates  the  currents  on  a sphere  of  "vacuum  di- 
electric . " In  view  of  the  equivalence  theorem,  the  cur- 
rents on  the  surface  should  obviously  be 
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One  notes  that  the  results  obtained  from  LDBR  code  agree 
well  with  the  expected  values.  Fig.  3.6  indicates  the 
currents  on  a dielectric  sphere  of  dielectric  constant 
er  * 2,  again  made  up  of  dummy  layers.  It  is  seen  that  the 
eigenfunction  solution  and  the  LDBR  solution  agree  favor- 
ably. The  results  are  also  in  good  agreement  with  those  of 
Wu  [14  ) and  Clisson  [ 11  ] , who  have  studied  the  homogeneous 
dielectric  sphere.  Fig.  3.7  indicates  the  equivalent  sur- 
face currents  obtained  from  LDBR  for  a homogeneous  sphere 
with  tr  = 1 and  various  values  for  the  conductivity.  The 
electric  currents  on  a perfect  electrically  conducting 
sphere  [l 5 ] are  plotted  for  comparison.  Fig.  3,8  indicates 
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Figure  3.8a.  Electric  currents  on  a sphere  of  high  dielectric  constant 
k a„  » 0.12,  k a_  - 0.2,  £ - 80,  O r*  0,  X = 1.0m. 
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the  currents  computed  using  LDBR  and  the  eigenfunction 
solution  procedure  for  a lossless  sphere  of  = 80,  to 
illustrate  the  capability  of  the  computer  code  to  treat 
"dummy"  layers  with  high  dielectric  constants.  Fig.  3.9 
and  Fig.  3.10  show  currents  computed  on  an  inhomogeneous 
sphere  with  different  medium  parameters  on  each  layer. 

The  results  of  LDBR  code  agree  very  well  with  those  obtained 
from  the  eigenfunction  solution  procedure.  It  should  be 
noted  that  the  values  of  conductivity  chosen  for  all  the 
cases  considered  here  are  rather  small.  This  is  because 
the  spherical  Bessel  functions  of  complex  argument  grow 
exponentially  with  increasing  conductivity.  In  fact,  with 
only  moderate  conductivities,  the  function  values  exceed 
the  dynamic  range  of  the  computer  used  to  perform  the  cal- 
culations, However,  the  code  LDBR  does  not  suffer  from 
this  limitation.  Although  the  LDBR  results  for  conductivities 
thus  cannotbe  checked  with  those  of  the  eigenfunction  formu- 
lation, nevertheless,  the  LDBR  results  do  converge  as  has 
been  found  by  increasing  the  number  of  sample  points  on  the 
layers.  Calculations  have  also  been  carried  out  for  the 
case  of  the  inhomogeneous  sphere  in  Fig.  3.9,  in  which 
additional  dummy  layers  have  been  inserted  within  the  homo- 
geneous regions.  It  has  been  found  that  treating  the  three 
region  problem  in  Fig.  3.9  as  a five-layered  sphere  (i.e. 


Figure  3.9b.  Magnetic  currents  on  the  inhomogeneous  sphere  of  Figure  3.9a. 


Figure  3.10a.  Electric  current  distribution  on  an  inhomogeneous 
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Figure  3.10b.  Magnetic  currents  on  the  inhomogeneous 
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two  layers  are  dummy  layers)  has  yielded  results  that  are 
in  excellent  agreement  with  those  obtained  using  but  three 
layers.  Fig.  3.11  shows  the  currents  on  a three  layered 
inhomogeneous  sphere  wherein  the  outermost  layer  is  chosen  to 
be  very  thin.  One  notes  that  the  results  of  the  LDBR  code  are 
still  in  good  agreement  with  those  of  the  eigenfunction 
solution  even  for  a body  made  up  of  thin  layers. 

We  finally  consider  some  additional  cases.  Fig.  3.12 
indicates  the  currents  on  a finite  homogeneous  cylinder  made 
up  of  dummy  layers.  Note  that  in  this  case,  the  dielectric 
region  has  surface  edges  and  hence  the  equivalent  surface 
currents  may  be  discontinuous  or  possibly  even  singular  at 
these  edges.  Fig.  3.13  indicates  the  currents  on  a homo- 
geneous dielectric  cone  sphere.  'Iso  shown  is  the  result 
obtained  by  Glisson  [16].  One  notes  immediately  the  differ- 
ence in  the  results  near  the  tip  of  the  cone.  This  difference 
is  attributed  to  the  layered  treatment  of  the  cone  sphere 
body,  which  introduces  a singularity  at  the  cone  tip  even 
for  source  and  field  points  not  on  the  same  surface.  Although 
the  current  expansion  scheme  we  adopt  does  not  match  the 
fields  at  the  tip,  nevertheless,  the  kernels  involved  will  be 
highly  peaked.  By  performing  the  Integration  around  the  cone 
tip  more  accurately,  this  error  is  easily  eliminated. 


istribution  on  a finite 


Figure  3.13a.  Electric  current  distribution  on  a 

dielectric  cone-sphere.  AAA  indicates 
results  obtained  with  improved  integration 
around  the  cone  tip. 
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CHAPTER  IV 

MISSILE  PLUME  SCATTERING 

In  this  chapter,  the  procedure  described  in  Chapter  III 
for  obtaining  the  scattered  fields  for  layered  dielectric 
bodies  of  revolution  is  extended  to  treat  a missile  with  an 
attached  inhomogeneous  exhaust  plume.  Calculated  missile 
currents  are  presented  for  selected  plume  inhomogeneities, 
f r equenc ies, and  angles  of  incidence. 

4 . 1 Block  Tridiagonal  Formulation  for  Missile  Plume  Scat- 
tering 

The  inhomogeneous  plume  is  modeled  as  a series  of  layers 
of  homogeneous  reg<ons.  Fig.  4.1  shows  an  approximate  model 
of  the  missile  plume  problem.  fine  boundary  conditions 
require  that  the  tangential  electric  and  magnoi  Lc  fields  be 
continuous  across  the  plume  layers  and  that  the  tangential 
electric  field  vanish  along  the  missile  surface  and  along 
the  missi le/ plume  interface.  This  leads  to  a set  of  coupled 
integral  equations  for  the  equivalent  currents  on  the  sur- 
faces between  layers  of  the  plume  and  on  the  missile  and 
mi ss i le / p lume  interface.  The  numerical  procedure 
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indicated  in  Chapter  III  leads  to  a matrix  equation  for  the 
representative  model  in  Fig.  4.1  as  follows: 
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to  the  field  contribution  on  surface  i due  to  currents  on 
surface  j.  An  immediate  observation  of  (4.1)  is  that  the 
block  tridiagonal  nature  of  the  overall  system  matrix  now 
appears  to  have  been  lost  with  the  presence  of  the  missile. 
However,  we  merely  regroup  some  of  the  operators,  currents 
and  driving  vectors  as  follows: 


L23  * ^L23  L24 
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32 

42 


L3  3 L34 

L43  L44 
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inc 

1 3 

V3 
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X4 
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_l5_ 

With  these  new  definitions,  (4.1)  may  be  written  as 
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which  retains  the  block  tridiagonal  nature.  The  extension 
to  an  arbitrary  number  of  layers  should  be  obvious.  Thus 
with  suitable  partitioning  of  the  matrix  equation,  the 
iteration  procedure  indicated  in  Section  2.4  can  still  be 
utilized. 

4 . 2 Boundary  Conditions  at  the  Plume-Missile  Junction 

The  boundary  condition  on  the  conducting  missile  is 
that  the  tangential  electric  field  vanish  over  the  missile 
surface.  On  the  plume  interfaces,  however,  we  require  that 
the  tangential  electric  and  magnetic  fields  be  continuous. 
Thus,  at  the  Junction  of  the  plume  and  the  missile,  we  must 
have  E * 0 and  to  be  continuous.  Since  J ■ n x H and 
M * E * n,  these  conditions  translate  into  conditions  on 
the  surface  currents  which  imply  that  the  t-directed 

component  of  the  electric  current.be  continuous  as  one 
approaches  the  junction  along  any  one  of  the  three  surfaces 
meeting  there  and  that  Mt>  the  t-directed  component  of  the 
magnetic  current,  vanish  at  the  junction  (Here  the  t- 
direction  is  the  direction  along  the  generator  arc  of  the 
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body  of  revolution).  Note,  however,  that  and  may  be 
singular  near  the  junction  since  these  currents  flow  paral- 
lel to  an  edge.  With  the  displaced  current  expansion  scheme 
(see  Fig.  4.2),  however,  these  current  components  are  not 
defined  directly  at  the  junction,  but  rather  a half- 
subdomain away  from  the  junction.  The  junction  point  of 
the  missile  and  plume  is  chosen  to  be  at  the  center  of  the 
t-directed  current  pulse  and  continuity  of  the  current  at 
the  junction  is  enforced.  Note  that  no  special  boundary 
conditions  on  or  at  the  missile-plume  junction  are 
needed  because  of  the  shifted  subdomain  scheme. 

In  order  to  apply  the  boundary  conditions  on  the  t- 
directed  component  of  the  electric  and  magnetic  currents, 
let  us  consider  the  equivalent  problems,  shown  In  Fig.  4.2, 
obtained  by  the  application  of  the  equivalence  principle. 

Fig.  4.2a  is  the  exterior  equivalence,  wherein  the  fields  in 
the  exterior  region  are  only  due  to  equivalent  electric  cur- 
rents residing  on  the  missile  surface  and  equivalent  electric 
and  magnetic  currents  residing  on  the  outside  of  the  plume 
surface.  The  fields  in  the  interior  are  assumed  to  be  zero 
for  this  part  of  the  problem.  Fig.  4.2b  shows  the  internal 
equivalence,  wherein  the  fields  in  the  exterior  region  are 
assumed  to  be  zero.  We  note  herein  that  the  fields  on  the 
interior  are  due  to  equivalent  electric  currents  residing  on 
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Figure  4.2b.  Internal  equivalence  and  current  expansion  sheme 
for  the  mi88ile/plume  configuration. 
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the  interface  of  the  missile/plume  and  due  to  electric  and 
magnetic  currents  on  the  interior  of  the  plume.  We  also 
note  that  the  missile  current  does  not  appear  in  this  part 
of  the  problem.  Using  the  above  two  equivalent  problems, 
we  note  that  the  following  are  true: 


<V  Jtj 

» “Ip  i 

Mp) 

+ h * Einc  « O' 

, (4.3) 

;(5M*  Jtj 

» Jp » 

Mp> 

+ n x fiinc  « 0, 

, (4.4) 

?£  S?  + SM 

“Etn<-V 

Jtj’ 

"V 

-Mp)  - O' 

) 

(4.5) 

Jtj’ 

-V 

-Mp)  - 0, 

> 

(4.6) 

re  Sp  + Sp 

where  Sp,  and  are  the  plume,  missile  and  interface 
surfaces,  respectively  and 


E is  the  electric  field  produced  by  the  currents 

*"A  for  the  exterior  equivalence,  evaluated  just 
inside  the  missile/plume  surface, 

H is  the  magnetic  field  produced  by  the  currents 

'"A  for  the  exterior  equivalence,  evaluated  just 
inside  the  miss ile /p lume  surface, 

E*n  is  the  electric  field  produced  by  the  currents 
for  the  interior  equivalence,  evaluated  just 
outside  the  p lume / int er f a ce  surface. 
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H^n  is  the  magnetic  field  produced  by  the  currents  for 
the  interior  equivalence,  evaluated  just  outside 
the  plume/interface  surface. 


is  the  electric  current  on  the  missile  surface. 


are  the  electric  and  magnetic  currents  on  the  plume, 


J , is  the  t-component  of  the  electric  current  at  the 
J junction, 


Jj  is  the  electric  current  on  the  missile/plume 
interface . 


Over  the  plume  surface  Sp,  (4.3)and  (4.5)  result  in 

a X r*  - fi  X E^n  - -ax  Einc  , re  Sp,  (4.7a) 


Similarly,  from  (4.4)  and  (4.6)  one  obtains 

a x H^x  - a X H^n  - -a  X iiinc  , reSp,  (4.7b) 

From  (4.3)  we  have, 

a X rx  - -a  X Einc  , reSM  , (4.8) 

From  (4.5)  we  have, 

6 X ®in  " 0 • ?eSi»  (4-9) 

where  r refers  to  the  point  at  which  the  fields  are  evaluated. 

Herein  the  dependence  of  the  various  fields  on  the  appropri- 
ate currents  are  suppressed.  Equations  (4.7)  through  (4.9) 
are  the  required  equations.  For  computational  purpose,  the 
layers  are  chosen  in  accordance  with  the  contours  of  the 
equivalent  problems  (see  Fig.  4.2).  If  one  uses  a testing 
procedure  similar  to  that  used  in  Chapter  III  for  the  layered 
dielectric  body  of  revolution,  the  above  equations  may  be 
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reduced  to  a matrix  equation.  At  this  stage  we  excercise 
some  caution  with  regard  to  the  equations  associated  with 
the  evaluation  of  the  fields  at  a point  due  to  various  cur- 
rent sub-domains.  For  test  points  on  the  plume  surface. 
Equation  (4.7)  essentially  corresponds  to  the  enforcement  of 
the  continuity  of  tangential  electric  and  magnetic  fields. 
Equations  (4.8)  and  (4.9)  correspond  to  the  vanishing  of  the 
tangential  electric  fields  over  the  missile  surface  and  inter- 
face , respectively.  At  the  junction,  the  net  effect  is  to 
add  together  three  equations  that  correspond  to  (i)  enforce- 
ment of  continuity  of  the  tangential  electric  field  at  the 
plume/exterior  region  interface,  (ii)  the  requirement  that 
the  tangential  electric  field  vanish  on  the  missile  surface, 
and  (iii)  the  requirement  that  the  tangential  electric  field 
vanish  at  the  missile/plume  interface.  The  use  of  the  cur- 
rent expansion  scheme  shown  in  Fig.  4.2  along  with  the  above 
testing  procudure  results  in  separate  terms  in  the  moment 
matrix  pertaining  to  the  t-directed  component  or  the  electric 
current  at  the  junction.  Since  these  terms  correspond  to 
the  same  unknown,  vi_z,  the  unknown  current  coefficient  at  the 
junction,  '-hey  are  grouped  together.  Such  a procedure 
yields  an  average  value  of  the  contribution  from  the  junc- 
tion current,  the  averaging  being  carried  out  over  all  the 
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three  surfaces  that  meet  at  the  junction  and  thus  achieves 
the  continuity  of  the  t-directed  component  of  the  electric 
current  at  the  junction.  Since  the  t-directed  magnetic 
current  pulses  overlap  the  t-directed  electric  current 
pulses  and  since  the  magnetic  current  on  the  missile/exterior 
region  interface  and  the  missile/plume  interface  is  zero, 
the  t-directed  magnetic  current  pulse  coefficient  at  the 
junction  is  forced  to  be  zero,  thus  yielding  a vanishing 
t-directed  component  of  magnetic  current  at  the  junc- 
tion . 

4 . 3 Numerical  Results 

The  computer  code  LDBR,  described  in  Chapter  III  for 
solving  currents  on  a layered  dielectric  body  of  revolution, 
vas  modified  to  incorporate  the  procudures  described  in 
Section  4.2  which  would  enable  one  to  solve  for  the  currents 
on  a composite  missile-plume  structure.  This  modified  code 
will  be  hereafter  referred  to  as  the  MPLM  code.  Since  the 
composite  missile-plume  body  is  tr.in  comparei'  to  the  length, 
only  the  circumferentially  uniform  Fourier  mode  (ra  = 0) 
current  is  considered.  This  is  because  for  reasonably  thin 
bodies  of  revolution  the  non-symmet r ic  modes  are  only  weakly 
coupled  to  the  plume  and  thus  these  modes  of  missile  current, 
though  perhaps  significant  by  themselves,  probably  do  not 
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vary  much  with  the  plume  presence. 

Due  to  the  lack  of  either  an  exact  solution  or  compar- 
able measured  results  for  the  missile/plume  geometry,  a few 
special  geometries  were  considered  and  the  results  compared 
with  those  obtained  from  other  valid  computer  codes.  Figs. 
4.3a  and  4.3b  show  the  currents  on  the  missile,  when  the 
plume  conductivity  is  set  to  zero.  As  one  would  expect 
upon  setting  the  plume  conductivity  to  zero,  the  currents 
induced  on  the  missile  should  be  identical  to  those  induced 
on  a conducting  cylinder  excited  by  a plane  wave.  Also 
shown  for  comparison  in  these  figures  are  the  results 
obtained  by  a computer  code  developed  by  Glisson  [16]  for 
the  currents  induced  on  a cylinder  due  to  a plane  wave  il- 
lumination. The  results  of  the  MPLM  code  are  in  excellent 
agreement  with  those  of  Glisson.  As  a further  check,  the 
geometry  of  a missile  with  a trailing  cylindrical  plume  (of 
the  same  radius  as  that  of  the  missile)  of  uniform  conduc- 
tivity 0 = 1000  S/m  is  considered.  In  this  case,  the  missile 
plume  combination  resembles  a perfectly  conducting  cylinder 
whose  length  is  the  sura  of  the  lengths  of  the  missile  and 
plume.  Figs.  4.4a  and  4.4b  show  the  circumferentially 
uniform  Fourier  mode  current  on  the  composite  structure. 

Also  shown  on  these  figures  are  the  results  obtained  by  the 
previously  mentioned  code  of  Glisson.  We  note  the  excellent 
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Figure  4.4b.  Electric  current  distribution  along  a raissil 
trailing  plume  of  uniform  conductivity 
O m 1000  S/m  under  oDlinut  incidence. 
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agreement  in  the  two  approaches.  The  excellent  agreement  in 
the  above  two  checks,  under  limiting  cases,  provide  some  con- 
fidence in  the  validity  of  the  computer  code  MPLM.  Calcula- 
tions were  also  made  for  a homogeneous  cylindrical  plume  with 
conductivity  0 - 0.2  S/m.  Fig.  4.5  depicts  the  resulting 
currents  on  the  missile-plume  combination.  Also  shown  are 
the  results  obtained  by  Wu  e t . a 1 . (17]  for  the  same  case. 

One  notes  a fairly  good  agreement  in  the  two  results.  A 
possible  explanation  for  the  difference  Is  the  different 
type  of  treatment  of  the  junction  between  the  missile  and  the 
plume  used  here  as  compared  to  that  used  by  Wu  e t . a 1 |17J. 
Further  it  is  known  that  the  junction  modeling  strongly 
Influences  the  currents  on  the  missile  (17],  Whereas  the 
procedure  n«»«d  here  essentially  averages  the  boundary  con- 
dition on  portions  of  the  missile,  plume  and  the  missile/ 
plume  interface  that  are  common  to  the  junction,  the 
approach  followed  by  Wu  e t ■ a 1 . [17]  is  t~  enforce  only  the 

boundary  condition  on  the  missile  and  let  the  enforced 
continuity  of  current  flowing  onto  the  teinainmg  surfaces 
at  the  junction  take  care  of  me  satisfaction  of  the 
boundary  conditions  on  these  surfaces.  This  latter 
approach  is,  in  principle,  correct;  however,  in  a num- 
erical procedure,  a certain  degree  of  "averaging"  of  the 
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boundary  conditions  is  desirable  because  a pulse  representa- 
tion of  the  current  is,  in  a sense,  a local  averaging  process 
and  hence  one  can  expect  the  boundary  conditions  to  be 
satisfied  only  in  an  average  sense  as  well.  Based  on  these 
observations  and  the  excellent  agreement  in  the  limiting 
cases  shown  in  Figs.  4.3  and  4.4,  it  is  concluded  that  the 
missile/plume  junction  treatment  used  here  probably  yields 
more  accurate  results  than  those  obtained  by  the  treatment 
adopted  by  Wu  et . al ■ {17],  and  that  the  differences  in  the 
results  seen  in  Fig.  4.5  are  due  to  the  different  ways  in 
which  the  junction  is  treated. 

For  the  actual  inhomogeneous  missile  plume,  the  plume 
inhomogeneity  is  predicted  using  the  LAPP  computer  code. 

This  code  simultaneously  models  thermo-chemical  reactions  and 
rocket  aerodynamics  to  establish  the  electrical  properties, 
viz.,  permittivity  and  conductivity,  inside  the  plume  region. 
The  plume  conductivity  calculations  we  have  used  here  model 
the  plume  of  a static  (zero  velocity)  Chapparal  missile  at 
a 5000  ft.  (1524m)  altitude.  The  calculations  show  that  the 
electrical  permittivity  does  not  change  much  from  the  free- 
space  value,  but  that  the  conductivity  is  strongly  inhomo- 
geneous both  radially  and  axially.  A detailed  discussion 
of  the  LAPP  code  is  given  in  [18].  Conductivity  values 


130 


predicted  by  the  LAPP  code  result  in  constant  conductivity 
profiles  shown  in  Fig.  4.6.  Based  on  this  accurate  profile, 
a judicious  choice  of  layer  boundaries  was  made  and  the 
assumed  constant  value  of  the  conductivity  between  the  layer 
boundaries  was  obtained  by  averaging  the  conductivity  between 
two  contours.  Fig.  4.7  shows  the  resulting  layered  approxi- 
mation to  the  inhomogeneous  plume.  We  shall  be  adopting  two 
models  of  the  plume,  the  long  and  short  plume  models.  Such 
a choice  is  made  for  two  reasons.  J.t  has  been  found  [17]  that 
the  effect  of  the  plume  on  the  currents  on  the  missile  is 
negligible  when  the  conductivity  in  the  region  around  the 
nozzle  of  the  missile  is  small.  Further,  the  values  of  the 
conductivity  in  the  plume,  as  predicted  by  the  LAPP  code,  is 
known  to  be  less  accurate  around  the  nozzle  region  as  compared 
to  the  values  predicted  in  the  regions  away  from  the  nozzle. 
The  extent  of  the  regions  of  the  short  and  long  plume  models 
are  shown  in  Fig.  4.7.  One  notes  that  the  region  around  the 
nozzle  has  a low  conductivity  value  in  the  long  plume  model 
and  a higher  conductivity  value  in  the  short  plume  model. 

Figs.  4.8  through  4.12  show  the  computed  currents  on  the 
missile  and  plume  for  various  angles  on  incidence  and  frequen- 
cies. Cost  of  computations  has  limited  a more  detailed  study 
of  the  short  plume  model.  We  note  from  these  figures  that  the 
currents  on  the  missile  are  not  much  affected  under  the  long 


Figure  4.7.  Layered  approximation  to  the  plume  of  the  Chaparral  missile. 


Figure  4.8c.  Electric  current  distribution  along  a Chaparral  missile  and  trailing 
plume,  f = 30  MHz,  0lnC  = 90°. 
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Electric  current  distribution  along  a Chaparral  missile  and  trailing 
plume,  f - 70  MHz,  einc  - 90°. 


d.  Electric  current  distribution  along  a Chaparral  missile  and  trailing 
plume,  f - 70  MHz,  0inC  - 120°. 


plume,  f - 90  MHz 


Figure  4. lie.  Electric  current  distribution  along  a Chaparral  missile  and  tra 


158 


plume  model  situation.  Due  to  the  low  conductivity  in  the 
plume  around  the  nozzle  region,  the  coupling  between  the  mis- 
sile and  the  plume  will  be  low.  Such  an  observation  has  been 
made  by  Wu  et . al . [17].  For  the  short  plume  model  case,  we 

note  that  that  the  nozzle  of  the  missile  is  now  in  a region 
of  a higher  conductivity  than  in  the  long  plume  model. The 
currents  on  the  missile  are  now  more  affected  by  the  presence 
of  the  plume.  One  notes  that  for  angles  of  incidence  grazing 
from  the  plume  side,  the  currents  on  the  missile  can  be  more 
than  the  no  plume  case.  One  notes  also  that  the  current  on 
the  missile  around  the  nozzle  region  is  much  higher  than  in 
the  no  plume  or  long  plume  cases.  Thus  coupling  into  the 
interior  of  the  missile  through  apertures  placed  in  the 
region  around  the  nozzle  of  the  missile  will  be  greater  in 
the  short  plume  case  than  in  the  long  plume  case. 

The  currents  on  the  plume  appears  to  be  comparable  to 
the  missile  current  at  low  frequencies.  However,  at  frequen- 
cies above  the  first  resonance,  the  current  on  the  plume  is 
generally  much  smaller  than  the  missile  current.  For  the  30MHz 
case  , (see  Figs.  4.8a  through  4.8e),  an  interesting  observa- 
tion in  the  plume  current  is  the  peak  occurring  at  about  3m 
from  the  junction  and  a shoulder  which  occurs  at  about  4.8m 
from  the  junction.  One  notes  that  the  skin  depths  at  this 
frequency  and  for  the  various  conductivities  involved  are  much 
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greater  than  the  radial  spread  of  the  plume.  Hence  assuming  a 
uniform  current  distribution  across  the  cross-section  of  the 
plume,  an  "equivalent  admittance",  Y , Per  unit  length  may 
be  computed  as  follows: 


Y 

eq 


> 


where  0 is  the  value  of  the  conductivity  in  the  elementary 
area  dA  at  the  cross-sectional  surface  S at  any  point  along 
the  axis  of  the  plume.  Fig.  4.13  shows  the  variation  of  the 
equivalent  admittance  along  the  plume.  One  notes  that  the 
location  of  the  peak  and  shoulder  in  the  equivalent  admittance 
variation  roughly  correspond  to  the  peak  and  shoulder  in  the 
current  on  the  plume. 
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Figure  4.13.  Equivalen 
of  the  pi 


t admittance  variation  along  the  axis 
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CHAPTER  V 

CONCLUSION 

In  this  work,  we  are  concerned  with  developing  an 
approach  for  treating  missiles  with  attached  inhomogeneous 
conducting  exhaust  plumes.  In  Chapter  II,  a discussion  of 
various  approaches  is  presented  and  a detailed  comparison  of 
the  unimoment  method  and  the  surface  integral  equation 
approach  with  a moment  method  solution  has  been  given.  It 
has  been  found  that  for  thin  and  layered  inhomogeneities, 
the  surface  integral  equation  approach  is  more  efficient, 
while  the  unimoment  method  is  more  efficient  for  highly 
varying  inhomogeneities  and  nearly  circular  objects.  Since 
the  objective  here  was  to  model  and  study  a composite  missile 
plume  configuration  wherein  the  axial  length  of  the  plume 
is  large  compared  to  its  radial  dimension,  the  surface  inte- 
gral equation  approach  was  adopted.  A computer  code  was  then 
developed  to  compute  currents  induced  on  layered  dielectric 
bodies  of  revolution.  The  validity  of  the  computer  code  was 
then  verified  by  comparing  the  computations  with  an  eigen- 
function formulation  for  concentric  layered  lossy  dielectric 
spheres.  The  surface  integral  equation  procedure  was  next 
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extended  to  the  missile-plume  problem  and  a simple  approach 
for  the  numerical  treatment  of  the  missile/plume  junction 
was  developed  and  validated. 

The  computations  on  the  layered  model  of  the  inhomo- 
geneous plume  have  yielded  results  which  are  in  good  agreement 
with  the  surface  impedance  boundary  condition  approach  used 
in  [17].  However,  the  surface  impedance  boundary  condition 
approach  should  be  viewed  with  caution.  This  is  because  the 
assumed  relationship  between  the  electric  and  magnetic  fields 
is  true  only  for  uniformly  illuminated  cylindrical  stuctures 
and  hence  can  be  expected  to  be  reasonably  valid  in  regions 
only  far  away  from  the  junction,  while  in  the  region  around 
the  junction  serious  errors  may  result  which  might  alter  the 
neighboring  current  distribution  along  the  missile.  Under 
such  circumstances,  the  layered  model  approach  should  yield 
more  accurate  results.  Cost  of  computations  has  precluded 
a more  detailed  analysis  of  the  short  plume  model.  Based  on 
the  good  correspondence  of  the  long  plume  model  results  ob- 
tained here  and  the  long  plume  model  results  obtained  in  [17], 
it  is  believed  that  the  currents  on  the  missile  under  short 
plume  model  conditions  would  be  more  strongly  affected  by  the 
presence  of  the  plume. 

In  this  work,  the  surface  integral  equation  formula- 
tions used  involve  both  the  electric  and  magnetic  field 
equations. 
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As  shown  in  Appendix  A,  these  problems  can  also  be  formulated 
in  terms  of  the  E-field  alone  or  H-field  equation  alone  or 
as  a combination  of  both  the  E and  H field  equations  to 
obtain  a so-called  combined  field  formulation.  These  ap- 
proaches have  been  applied  here  only  to  homogeneous  dielec- 
tric cylinders.  However,  they  may  be  extended  easily  to 
three-dimensional  bodies  and  may  provide  useful  alternative 
approaches . 

One  worthwhile  approach,  which  has  not  been  considered 
here  in  the  evaluation  of  methods  for  treating  the  missile/ 
plume  problem,  is  to  couple  an  integral  equation  on  the 
missile  and  plume  boundaries  with  a finite  element  or  finite 
difference  solution  of  the  fields  in  the  inhomogeneous  plume 
region.  Such  an  approach,  which  has  so  far  not  been  carried 
out,  combines  the  best  aspects  of  both  the  unimoment  method 
and  the  surface  integral  equation  approach  and  should  be 
applicable  to  rather  arbitrary  shaped  geometries  and  inhomo- 
geneities. Furthermore,  the  approach  also  appears  as  an 
attractive  method  for  analyzing  coupling  through  apertures 
into  the  interior  of  missiles  or  other  structures.  In  :uch 
problems,  the  exterior  region  is  formulated  by  a surface 
integral  equation,  while  the  interior  is  formulated  in  terms 
of  the  appropriate  wave  equation.  The  two  equations  are 
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then  coupled  at  the  aperture  with  the  appropriate  boundary 
conditions . 
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APPENDIX  A 

SURFACE  CURRENT  FORMULATIONS  FOR 
DIELECTRIC  SCATTERING  PROBLEMS 


Scattering  by  homogeneous  dielectric  bodies  can  be 
formulated  via  the  surface  equivalence  principle.  Applica- 
tion of  th$  boundary  conditions  leads  to  a set  of  four 
integral  equations  involving  the  two  unknown  equivalent 
surface  currents  J and  M.  However,  under  suitable  conditions, 
only  two  of  the  equations  are  sufficient  to  determine  J and 
M.  Thus  various  combinations  of  these  four  equations 
can  be  used,  each  combination  leading  to  a different  type 
of  formulation.  Some  of  the  more  important  choices  are 
described  by  Mautz  and  Harrington  [10]  . One  common  choice 
of  the  combination  leads  to  the  so-called  PMCHW  formulation 
described  in  [10].  Mautz  and  Harrington  also  point  out  that 
the  so-called  Mueller  formulation  is  related  to  the  PMCHW 
formulation  in  that  both  approaches  are  special  cases  of  a 
"combined  field  formulation."  We  consider  here  those  com- 
binations and  further  point  out  some  features  of  combina- 
tions other  than  those  considered  in  [10].  For  the  sake  of 
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clarity  and  completeness , some  of  the  material  in  [10]  is 
repeated  here. 

a . 1 Surface  Integral  Equation  Formulation 

Fig.  A.l  depicts  the  homogeneous  scatterer  illuminated 

i n c ••inc 

by  an  incident  field  (E  , H ).  The  body  parameters  are 
denoted  by  (yd,  ed> , the  fields  inside,  by  (Ed>  Hd> . The 
body  is  immersed  in  a medium  characterized  by  parameters 
(p^,  ee).  The  straightforward  application  of  the  equivalence 
principle  leads  to  the  exterior  and  interior  equivalences 
shown  in  Figs.  A. 2 and  A. 3 respectively.  Equating  to  zero 
the  tangential  components  of  the  null  fields  appearing  in 
Figs.  A. 2 and  A. 3 leads  to  the  following  equations: 


-n 

x Eg(J,  M)  » n * 

EinC 

> 

(A.l) 

-A 

x H"U,  M)  - n x 
e 

HlnC 

» 

(A. 2) 

-n  * Ed ( J , h)  - 

0 , 

(A. 3) 

-n  x Hd(3,  M)  - 

0, 

(A. 4) 

where 

V “i  are 

the  electric  and 

magnetic 

fields  due  to  J 

and  M 

radiating 

in  a medium 

e^)  and  + 

indicates 

Figure  A. 1.  A homogeneous  dielectric  scatterer. 


a 


s 


Figure  A. 3.  Internal  equivalence. 
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evaluation  of  the  fields  Just  inside  and  outside  S,  respec- 
tively. We  note  that  (A.l)  - (A. 4)  is  a set  of  four  equa- 
tions in  the  two  unknowns  J and  M.  We  need  to  reduce  the 
number  of  equations  to  two  and  this  can  be  done  by  ap- 
propriately selecting  or  combining  the  four  equations  in 
various  ways.  Some  of  the  various  possible  combinations  of 


(A.l)  - (A. 4)  are 

as 

follows i 

(a)  -fi  x 

+ alt)  ■ ft  x 

d 

sine 
b » 

(A. 5) 

-ft  x 

<s; 

+ BH*)  -fix 

HlnC  , 

(A. 6) 

(b)  -fi  x 

h" 

p 

-fix  HinC  , 

(A. 7) 

■fi  x 

a 


(c)  -ft  x E 


0 , 


ft  x EinC 


•n  x Ej  « 0 

a 


(A. 8) 


(A. 9) 


(A. 10) 


(d)  -ft  x fi"  - — E“  . - ft  x HinC  + ~ E*nC  (A. 11) 

e ng  e,tan  He  tan 


n x fi*  + it 


d nd  d , t a n 


0 , 


(A. 12) 
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where  ot  and  8 are  arbitrary  complex  constants.  The  first 
type  of  combination  in  (a)  above  has  been  considered  by 
Mautz  and  Harrington  [10].  As  they  have  pointed  out,  the 
choice  of  a * 8 * 1 leads  to  the  PMCHW  formulation,  while 
the  choice  of  8 * “ (Vj/Vg)  leads  to  the  Mueller 

formulation.  We  designate  combination  (b)  as  the  H-field 
formulation  (HFIE) , combination  (c),  the  F*field  formulation 
(EFIE),  while  (d)  is  called  the  combined  field  formulation 
(CFIE) . It  should  be  noted  that  the  CFIE  formulation  pre- 
sented here  in  (d)  is  different  from  the  formulation  (a) 
above,  which  was  also  called  a ''combined  field  formulation" 
by  Mautz  and  Harrington  [10].  Our  CFIE  formulation  (d)  is 
actually  a generalization  of  their  so-called  combined  field 
formulation  for  perfectly  conducting  scatterers  [19].  Note 
that  (d)  combines  two  types  of  field  quantities  in  each 
equality,  whereas  (a)  combines  the  same  type  of  fields  from 
different  regions.  We  consider  in  detail  the  combinations 
(b)  - (d). 

A.  2 H-Field  Formulation 

The  basic  equations  for  this  type  of  formulation  are 
-n  x H~( J , M)  = n * Hlnc 
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-A  X H*(J,  M)  - 0. 


Here,  we  have 


H^j,  M)  * ~ j toF1  (M)  - V<f>“(M)  + ~ V x A^J),  (A. 13) 

i ■ d or  e , 


where  F^M)  is  the  electric  vector  potential  due  to  the 
magnetic  current  M in  a medium  characterized  by  the  para- 
meters (y^,  $”(M)  i 8 the  magnetic  scalar  potential  due 

to  M in  a medium  characterized  by  the  parameters  (y^,  E ^ ) 
and  Ai(J)  is  the  magnetic  vector  potential  due  to  the  elec- 
tric current  in  a medium  characterized  by  the  parameters 
(*Ji,Ei>*  ExPresaing  the  ootentials  in  terms  of  the  integrals 
over  J and  M,  one  obtains 


J 


J 


+ L 


HH 


+ L 


HH 


(M)  - p x HlnC  , (A. 14) 

(M)  - 0 f (A. 15) 


where 
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•ffj  x VG  t 
S 


ds 


-Jwei 


jj  MGids' 
S 


V 


(V'  • M) 
(- jw) 


ds' 


i =■  d,e. 


where  G^  is  the  appropriate  homogeneous  Green's  function 
for  the  medium  with  the  parameters  ( ^ The  notation 
for  the  operators  above  is  as  follows:  the  first  subscript 
refers  to  the  evaluation  of  the  type  of  field,  which  in  this 
case  is  the  H- field,  while  the  second  subscript  refers  to 
the  type  of  source,  either  electric  (indicated  by  E) , or 
magnetic  (indicated  by  H) . The  superscripts  d and  e refer 
to  the  medium  characterized  by  (p^,  C^)  and  (Ue*  ee)> 
respectively.  The  factor  1/2,  where  I denotes  the  identity 
operator,  results  from  the  evaluation  of  V x A on  the  sur- 
face of  the  scatterer.  One  may  now  use  (A. 14)  and  (.‘  IS) 
to  solve  for  J and  M via  moment  method.  An  alternative  ap- 
proach would  be  to  use  (A. 15)  to  express  M ia  terms  of  J as 

B ■ <lhh)_1  <1  - 4>'J 


(A. 16) 
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Using  (A. 16)  in  (A. 14),  one  then  obtains 

[l  + lHE  + lHH  <1  - I4e)]  3 * * * S1""'  <A'17> 

Equation  (A. 17)  can  be  used  to  solve  for  J,  Having  deter- 
mined J,  M can  then  be  obtained  from  (A. 16).  We  point  out 
that  by  storing  the  matrix  (L^)  1 (-|  - L^)  , which  is 
computed  in  the  process  of  determining  J,  one  may  compute 
M directly  once  J has  been  determined. 

It  should  be  noted  that  expressing  M in  terms  of  J, 
as  in  (A. 16),  presupposes  the  existence  of  (L^)  How- 
ever, it  can  be  shown  (see  (10])  that  does  not  possess 

rln 

an  inverse  at  frequencies  corresponding  to  resonant  fre- 
quencies of  a conducting  cavity  formed  by  the  closed  surface 
S and  filled  with  a material  having  electrical  parameters 
(p.,  e,).  One  may,  however,  solve  for  J and  M from  (A. 14) 
and  (A. 15)  simultan-  Dusly  at  such  frequencies  (10],  We 
introduce  here  a terminology  to  describe  these  two  ap- 
proaches. Whenever  J and  M are  simultaneously  solved  as  in 
(A. 14)  and  (A. 15),  we  call  the  approach  the  two  current 
formulation,  while  the  use  of  (A. 16)  and  (A. 17),  in  which 
one  solves  for  one  current  at  a time,  is  called  the  single 


current  formulation. 
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A.  3 E-Field  Formulation 


The  basic  equations  in  this  approach  are 


a v p"  _ a =inc 

“ n * — n * b 

e 


-n  x E*  - 0, 


Here  we  have 


V4>?( J)  - ~ V X F (fi) 
1 i 1 


i » d or  e . 


(A. 18) 


Expressing  the  potentials  in  terms  of  the  currents,  we  have 


lee(3)  + 1 + leh  fi  ' a - ilnc 


lee<3)  + leh  - 1 R ■ 0 ■ 


(A. 19) 


(A. 20) 


where 


lee(j) 


tt  - rrSTill  !a 

•jw^i  jl  JG^ds ’ - Vjj  <-ju>)  (Ji 


ds’  , 
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LgR(M)  * x 7Gj[  ds’  , 

S 

with  the  notations  for  the  operators  paralleling  that  of  the 
H-field  formulation.  Proceeding  in  a manner  similar  to  the 
H-field  formulation,  one  obtains  the  single  current  E-field 
formulation  as 


0 + '■EH  + 4(14> 


-1 


- L 


eh]]  fi 


x £inC 


(A. 21) 


with 


3 ■ (LEErl  (I  - 4 )S  • 


(A. 22) 


As  in  the  case  of  the  H-field  formulation,  the  single  current 
E-field  formulation  fails  at  the  interior  resonant  cavity 
frequencies.  However,  as  before,  one  may  adopt  the  two 
current  approach  of  solving  for  J and  M s imultaneously  from 
(A. 19)  and  (A. 20)  at  such  frequencies. 


A . 4 Combined  Field  Formulation 

The  basic  equations  in  this  formulation  are 


H - 
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'e , tan 


n x 


HinC  + 


JL  FinC 

n tan  ’ 
e 


175 


' 


SO 


-fi  x iij  + 1 e* 

d rij  d,tan 


where  r)  and  n.  are  the  characteristic  impedances  of  the 
e d 

exterior  and  interior  medium  respectively.  Expressing  the 
fields  in  terms  of  the  currents  J and  M,  we  obtain  the 
following  equations; 


0,  (A. 24) 


where  the  operators  are  defined  as  before.  A feature  of 
this  formulation  is  that  the  fields  have  been  combined  in 
a manner  identical  to  that  of  the  combined  field  formulation 
for  conducting  bodies  proposed  by  Mautz  and  Harrington  [19], 
This  fact  makes  it  relatively  easy  to  incorporate  a solution 
procedure  for  both  a dielectric  object  and  a conducting 
object  within  the  same  computer  program.  This  can  be  easily 
seen  by  noting  that  (A. 23)  becomes  the  combined  field  inte- 
gral equation  for  a conducting  scatterer  [19],  if  one  simply 
retains  only  terms  involving  the  electric  current  on  the 


3 


i c 

i t 


left  hand  side  of  the  equation.  Following  [10],  the 
permissible  values  of  a and  8 can  be  determined  by  con- 
sidering for  what  values  of  a and  8 the  homogeneous  equa- 
tions 


-n  x H-  - — E~  . * 0, 

e n e , tan  * 
e 


(A. 25) 


A -+  8 — I 

-n  x H,  + E 


d nd  d,tan 


0, 


(A. 26) 


have  only  the  trivial  solution.  The  first  step  in  doing 


this  is  to  show  that  E . , H~  _ , Et  and  H+ 

e.tan’  e.tan  d,tan  d,tan 


are  all  zero.  For  this,  we  consider  first  the  complex 
power  flow  into  the  interior  of  surface  S in  the  equivalent 
problem  (A. 25) : 


Pe  * ff*e  * fie  * (_fi)  ds 


JJr  . <n  x rr*)  d8 


a 

n 


// 


E ds 

e tan  1 


where  (A. 25)  has  been  used  in  (A. 27).  We  consider  two  cases 


(1)  The  medium  with  parameters  yg,  eg  is  lossless. 

Then  Re  (P  ) = 0.  But  from  (A. 27),  if  Re(— ) 4 0, 

ne 


we  note  E = H 

e,tan  e,tan 


0. 


(2)  The  medium  with  parameters  y , e is  lossy.  If 

6 G 

Re(Pg)  * 0,  then  by  the  uniqueness  theorem  [13], 

E * H * 0 interior  to  S.  If  we  assume  Re(P  ) > 0 
e e e 


we  conclude  E 


" - H~  „ - 0.  If  Re (—-)  - 0, 

e,tan  e,tan  n 


then  Re(P  ) =*  0 , and,  as  above,  we  obtain  E 

e e , tan 

- H~  „ - 0. 

e , tan 


The  argument  above  is  also  valid  for  the  perfect  conductor 


case  and  extends  the  allowable  values  of 


beyond  those 


cos idered  in  [19]. 

Next  consider  the  complex  power  flow  into  the  exterior 


of  S in  the  equivalent  problem  (A. 26): 
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• n ds 


(n  x H*)*ds 


tan 


2 


ds 


(A. 28) 


where  (A. 26)  has  been  used  to  obtain  (A. 28). 


If  Re(p  ) - 0,  then  E, 
a a 


■ 0 by  the  uniqueness 


theorem  [13).  Thus,  assume  Re(P^)  > 0.  Then  if  we  choose 


R e 


8 


> 0,  we  conclude  E 


6 


. ■ H , =0  from  (A. 28) . 

d , tan  d , tan 


'd' 


“ 0,  we  have  Pe(P  ,)  » 0 and  from  the 

a 


H 


0. 


Further,  if  Re 
above  argument  E^  - 

To  show  next  that  E 

e,  tan 

implies  J = M * 0,  we  consider  the  exterior  aquivalence  shown 


Ej  **  H ■ Hj 

d,tan  e,tan  d,tan 


in  Fig.  A. 4,  wherein  E 


e , tan 


H * 0,  but  where  we  assume 

e , tan 


E£,  H outside  the  scatterer  to  be  non-zero.  Similarly  for 
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the  interior  equivalence  shown  in  Fig.  A. 5,  we 

* H+  _ « 0,  but  we  assume  E,,  ii,  interior  to 

to  be  non-zero.  Since  the  interior  fields  in  Fig.  A-4  are 
zero,  we  may  change  the  internal  medium  to  (p^ , e^) (see  Fig. 
A. 6),  Similarly,  in  Fig.  A. 5,  the  exterior  fields  are  zero 
so  that  we  may  change  the  exterior  medium  to  (Me»£e)  (See 
Fig.  A. 7).  We  note,  however,  that  Figs.  A, 6 and  A. 7 depict 
identical  situations  with  regard  to  sources  and  media  except 
for  the  change  in  the  sign  of  the  sources.  Hence,  except 
for  the  sign,  the  fields  radiated  by  the  currents  should  be 
the  same  in  both  figures.  Hence  (Ee,He>  = (Ed>H^)  * (0,0), 
which  implies,  in  turn,  that 

J - n x H+  - 0 
e 

M - _ n x E+=  0 
e 

Thus  we  have  shown  that  (A. 25)  and  (A. 26)  possess  only 
trivial  solutions  and  hence  that  the  solutions  of  (A. 23)  and 
(A. 24)  are  unique . 

One  may  present  a circuit  analogy  which  helps  to  explain 

the  conditions  on  a and  8 . Referring  to  Figs.  A . 8 =nd  A . 9 ^ 

we  may  think  of  Ye  - — and  Y1^  = as  surface  impedances 

s n s n, 

e d 


have  it  _ 

d , tan 


the  scatterer 


Figure*  A.  4.  Exterior  equivalence  I’lgure  A.$.  Interior  equivalen  s 
with  assumed  non-ssero  with  asHumcd  non-."  ro 

fields  In  the  exterior  fields  in  the 

region.  interior  region. 


Figure  A. 6.  Exterior  equivalence 
with  interior  medium 
parameters  (y^,  e^). 


Figure  A. 7.  Interior  equivalence 
with  exterior  medium 

parameters  (y  , e }. 

e r 
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and  (A. 25)  and  (A. 26)  as  surface  Impedance  boundary  condi- 
tions. The  problem  of  determining  whether  or  not  there 
exist  non-trivial  solutions  to  (A. 25)  and  (A. 25)  is  thus 
one  of  deciding  whether  or  not  the  internal  or  external 
regions  can  support  a resonance  when  the  region  is  bounded 

g d 

by  the  surface  impedance  Y or  Y , respectively.  For 

S 8 

example,  in  the  external  resonance  problem,  since  there  is 
always  some  loss  due  to  radiation,  then  if  the  surface  is 

g 

reactive  or  has  a small  loss  (Re(Y  ) _>  0)  , there  can  be  no 
resonances.  For  the  internal  problem,  if  the  medium  is 

g 

lossless,  then  no  resonances  are  possible  if  Re(Y  ) 4 0, 

8 

g 

since  the  surface  is  either  lossy  (Re(Y  ) > 0)  and  therefore 

s 

the  fields  are  damped,  or  supplies  power  to  the  Interior 

region  (Re(Ye)  < 0)  and  therefore  the  fields  grow  with  time 
s 

since  there  is  no  corresponding  absorpcion  mechanism.  In 
the  lossy  case,  however,  we  may  not  have  internal  resonances 

g 

if  Re (Y  ) > 0 because  of  internal  and  surface  losses, 
s — 

A. 5 Application  of  Various  Surface  Integral  Formulations 
to  TM  Scattering  by  Dielectric  Cylinders 
To  illustrate  some  of  the  foregoing  observations,  we 
use  the  method  of  moments  with  some  of  the  various  formula- 
tions discussed  here  to  solve  scattering  from  a homogeneous 
dielectric  cylinder  excited  in  the  TH  polarization.  Assuming 


war 


Figure  A. 8.  Circuit  analogy  and  permissible  values  of  a/ne 
for  exterior  equivalence. 


Permissable 
values  of  e/n. 


in(AA^) 


9»»  r*(/Ja^) 


Figure  A. 9.  Circuit  analogy  and  permissible  values  of  3/n^ 
for  interior  equivalence. 
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that  the  axis  of  the  cylinder  is  along  the  z-directlon  (Fig. 
A.  10) , we  have 


3 « J £ , (A. 29) 

z 

M - Mx  t , (A. 30) 

We  note  also  that  V'  • J ■ ,0,.  since  ,the  cylinder  is  Infinite 
in  the  z-direction.  The  various  operators  can  now  be  writ- 
ten as 


LEE<3)  * t / ViHo2)  <**  Ip-P'I)  dc'  , (A. 31) 


leh(H>  • Tj 


/ \ Ti’  "o2>(kl  l°-5' I)  dc',  (A-32) 


4e<3)  ' Tf  i i IP-P’I)  dc'.  (A-33) 


e;h(M)  - j / dc' 


_1_  _9_ 

4u  9 x J aT 
c 


/ 


?Mx  K^^kjp-p‘1) 


dc ' , 


(A. 34) 


-***fS3^  V^JT 


Figure  A. 10.  Geometry  of 
cylindrical 
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1 / 2 ) _ ^ 

where  the  two  dimensional  Green's  function  tF  y ( k ^ | p — p * | ) 
for  the  medium  with  parameters  C y ^ ^ ) has  been  used,  and 


,(2) 


k,  » (j/pTeT,  i «■  d or  e.  Hv“'(x)  is  the  Hankel's  function 
l i i o 

of  second  kind  and  zero  order,  p and  p'correspond  to  field 
and  source  points,  respectively. 

In  order  to  numerical ly  approximate  these  operators, 
we  divide  the  contour  of  the  cylinder  into  a number  of 
straight-line  segments  as  shown  in  Fig.  A. 11.  Pulse  func- 
tion basis  sets  are  used  to  represent  each  of  the  currents 
over  the  subdomains.  Thus  the  current  is  expanded  as 


N 


J (t)  ■ J p (t) 

z ^ n*n 


n - 1 


(A. 35a) 


N 


Mt(£)  ■ 2 Vn(t)  • 


n ■ 1 


(.A.  35b) 


where 


pn(o 


‘1,  t , < t < t , . 

’ n-^  — — n+*5 


.0,  otherwise 


(A. 36) 


and  t is  the  arc  length  along  the  cylinder  contour.  Using 
(A. 35),  the  vector  potential  contributions  in  (A. 31)  - (A. 34) 


im-i 
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can  be  computed.  In  order  to  compute  the  scalar  potential, 
we  use  a scheme  similar  to  that  of  Glisson  [11],  wherein  a 
pulse  representation  of  the  charge  is  derived  from  a finite 
difference  approximation  of  the  continuity  equation.  Thus 
we  have 


N + 1 


j 3H 

p.(tl  ■ i 17 


42 


M - M 


n-1 


n-*s 


(t) 


(A. 37) 


n-1 


where  the  charge  pulses  are  defined  as 


pn-* 


(t) 


1. 

0, 


t , < t < t 
n-1  — — 

otherwise 


n 


(A. 38) 


A point  matching  testing  procedure  is  used  to  evaluate  the 
operators  at  the  center  of  current  pulses.  The  gradient 
term  in  (A. 34)  is  evaluated  by  a finite  differencing  pro- 
cedure as  in  [11).  By  defining  the  "total"  arc  points,  t^ 
of  the  contour  to  be  at  possible  bends,  the  fields  wilx 
always  be  matched  away  from  points  wherein  fields  may  be 
singular.  With  the  above  expansion  and  testing  procedures, 
the  various  operators  may  be  approximated  as  follows: 


where 


iSf-jw;?51  ffr-  i~'_w  wn 
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'3  4wy 


n 

hr  I Ho2)(tiiw;i) 

oun  * 


n-1 

t 


n+1 


'4  4wy 


±T S H^2>  Ck±[p„  v-p;  I > <!<=• 

o4n+l  * ° 1 “ * " 


C,  - - 


5 4wy . 


(!  - 1 i l»  -y'  * dnr  <*«  t - u 


vV 

I *j 


'6  4u,y 


(1  ■ j ! *n  ^)  + t-*—  (£n 


n+1 


it  2 2irwy  v~“  2 

o 


- 1). 


Computer  program  subroutines  for  the  evaluation  of  each 
of  the  operators  Lgg,  Lgy,  bgE  «nd  LHH  were  written  and 
combined  in  the  appropriate  manner  depending  on  the  type  of 
formulation  used.  In  Fig.  A. 12  is  shown  the  currents  on  a 
circular  cylinder  as  found  by  the  various  methods.  Also 
shown  is  the  exact  eigenfunction  solution.  Fig.  A. 13 
depicts  the  currents  on  a square  cylinder  of  side  2a.  Since 
the  single  current  formulation  leads  to  erroneous  results 
at  the  resonant  frequencies  of  the  interior  region,  a plot 
of  the  determinant  of  the  moment  matrix  vs.  ka  is  shown  in 
Fig.  A. 14.  We  note  that  both  E and  H types  of  formulations 
indicate  a sharp  dip  in  the  value  of  the  determinant  at  the 
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resonant  frequencies,  which  are  the  cut  off  frequencies  for 
modes  of  the  square  waveguide  filled  with  a dielectric  with 
er  * 2.56.  The  combined  field  formulation  when  solved  by 
the  single  current  approach,  also  shows  the  resonance  effect. 
This  Is  because  the  operators  on  J and  M,  which  otherwise 
could  be  used  to  express  J in  terms  of  M or  vice  versa,  do 
not  possess  an  inverse  at  these  resonant  frequencies.  How- 
ever, when  one  solves  (A. 23)  and  (A. 24)  simultaneously  one 
avoids  this  difficulty  as  shown  by  the  plot  of  the  deter- 
minant of  the  CFIE  matrix. 
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